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1  Introduction 


Slater  orbitals  provide  a  well-known  basis  for  efficient  descriptions  of  the  atomic  and  molecular 
eigenstates  required  in  studies  of  chemical  structures  and  related  properties,  and  are  potentially 
applicable  to  a  wide  range  of  AFRL  and  other  DoD  branch  interests.  A  recent  series  of  studies  [HI 
111  describes  a  computational  code  suite  (SMILES)  developed  in  our  group  for  calculations  of 
molecular  integrals  which  overcomes  many  of  the  difficulties  that  have  previously  precluded 
wide-spread  adoption  of  Slater-based  methodology  for  these  purposes. 

The  SMILES  methodology  in  present  form  is  largely  applicable  to  the  orbital  principal  quantum 
numbers  required  for  accurate  descriptions  of  the  valence  shells  of  atoms  and  molecules,  and  is 
therefore  highly  suitable  largely  for  studies  of  ground- state  molecular  structure  and  low-lying 
valence  excited  electronic  states.  Highly  excited  molecular  electronic  states,  however,  are  of 
considerable  interest  in  a  great  many  connections,  including  optical  spectroscopic  diagnostics 
for  identification  and  characterization  of  newly  synthesized  chemical  compounds  for  propulsion 
and  energetics [O,  studies  of  photochemical  potential  energy  surfaces  for  atmospheric  and  re¬ 
lated  reactions 01,  descriptions  of  molecular  Rydberg  and  continuum  states  required  in  photon 
and  electron-impact  excitation  processes  0,  and  in  calculations  which  require  spectral  closure 
over  electronic  states,  such  as  the  long-range  interactions  required  in  atomic  and  molecular 
cluster  studies  fl6j|,  to  mention  some  representative  examples. 

Extension  of  SMILES  methodology  to  higher  principal  quantum  numbers  would  open  the  way 
to  development  of  a  great  many  new  computational  applications  suites  which  can  provide  con¬ 
comitant  support  to  on-going  DoD  materials,  chemical  physics,  and  aeronautical  and  propulsion 
sciences  research  programs.  This  extension  requires  a  thorough  revision  on  the  performance  of 
the  algorithms  currently  available  and  the  design  of  new  alternatives  for  those  cases  in  which 
the  current  algoritms  are  not  well  suited  for  large  quantum  numbers. 

The  currently  available  SMILES  methodology  is  satisfactory  for  atomic  orbital  involving  usual 
values  of  screening  factors  and  angular  momentum  quantum  numbers,  L,  ranging  from  0  to  5, 
with  allowable  principal  quantum  numbers,  N,  from  1  to  7  —  L. 

The  methodology  has  been  tested  in  a  high  amount  of  molecular  calculations  dealing  with  the 
ground  and  low-lying  excited  states  in  molecules  for  standard  equilibrium  geometries.  However, 
the  performance  of  these  algorithms  in  less  usual  situations,  as  those  requiring  very  high  or  very 
small  screening  constants  and  higher  quantum  numbers,  has  not  been  examined  to  date. 

The  current  work  starts  with  the  thorough  analysis  of  the  performance  of  the  algorithms  cur¬ 
rently  implemented  in  SMILES  in  a  wide  range  of  values  of  screening  constants  and  for  mod¬ 
erately  high  values  of  quantum  numbers.  In  those  cases  in  which  the  algorithms  do  not  provide 
sufficient  accuracy,  other  alternatives  are  formulated  and  coded,  and  their  performance  is  tested 
in  order  to  assess  their  capability  to  overcome  the  existing  ones. 

The  work  has  been  developed  in  steps  corresponding  to  the  cathegories  of  integrals  appearing 
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in  molecular  calculations.  In  each  step,  a  robust  code,  usually  lengthy  for  practical  pruposes  but 
good  enough  for  a  reference,  has  been  developed  to  provide  a  reliable  reference  in  testing. 
Once  the  code  for  the  reference  has  been  implemented,  integrals  corresponding  to  appropriate 
quantum  numbers  have  been  computed  for  a  wide  set  of  screening  parameters.  The  currently 
available  algorithms  and  the  reference  program  have  been  used  for  this  purpose,  and  a  thorough 
comparison  on  accuracy  has  been  carried  out.  In  those  cases  in  which  the  current  algorithms 
have  revealed  unsatisfactory,  alternative  procedures  have  been  formulated,  coded  and  tested. 

In  all  cases,  even  for  the  algorithms  already  available  in  SMILES,  new  codes  in  FORTRAN  90 
have  been  developed  in  order  to  allow  their  extension  to  quadruple  precision  and  multiprecision 
(mainly  for  the  reference).  In  some  cases,  this  has  implied  drastic  modifications  in  the  codes 
with  respect  to  the  original  ones,  with  new  algorithms  for  the  computation  of  auxiliary  functions 
and  a  significant  storage  reorganization. 

In  the  following  sections,  after  a  summary  of  the  basic  definitions  and  concepts,  the  fundamen¬ 
tals  of  the  techniques  applied  for  solving  the  integrals  are  described  and  the  algorithms  analyzed 
for  each  type  of  integrals  are  reported.  To  make  the  text  more  readable,  results  and  conclusions 
are  presented  and  commented  within  each  type  of  integrals. 

The  report  ends  with  a  list  of  the  programs  developed.  The  files  containing  the  source  codes  can 
be  found  as  supplementary  material  accompanying  this  report. 


2  Basic  concepts 


LCAO  ab  initio  calculations  of  the  electron  structure  of  molecules  require  four  types  of  inte¬ 
grals:  overlap,  Srs,  kinetic  energy,  Trs,  nuclear  attraction,  V/s,  and  electron  repulsion,  vrstu. 

For  a  given  basis  set:  { Ar- } "L  ] ,  they  are  defined  respectively  by: 


Srs  =  dr  Xr( r)  Xs(r)  =  (  Xr  |  X«  ) 


T  =  — - 
rs  2 

V*  =  -Qi 


dr 


Xr(r)  Xs(r)  = 
|r  —  R/| 


where  Qj  is  the  charge  of  nucleus  I,  and 

Vrstu  =[  dr  Xr{ r)  X«(r)  [  dr' 


r  —  r 


!  Xr  I  Xs  ) 

(1) 

=  (Xr\T\Xs) 

(2) 

Xr  |  |  Xs  ) 

rj 

(3) 

,(r0  r  |  n 

h  —  [  Xr  Xs  \  Xt  Xu  \ 

(4) 

As  the  definitions  render  evident,  molecular  integrals  can  be  classified  into  two  main  cathegories 
attending  the  dimension  of  the  integrals:  one-electron  (3D)  and  two-electron  (6D).  It  is  also  well 
known  that  the  difficulties  found  in  the  evaluation  of  integrals  with  STO  strongly  depend  on  the 
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number  of  centers  of  the  functions.  According  to  this  number,  they  are  further  divided  in  one, 
two,  three  and  four-center  integrals. 

This  report  deals  with  the  calculation  of  molecular  integrals  involving  exponential  functions 
known  as  Slater  type  orbitals  (STO).  An  unnormalized  STO  centered  at  a  point  R,i  is  a  product 
of  a  regular  (real)  harmonic  times  a  radial  factor: 

(5) 

where  =  r  —  R  4,  n  is  a  nonnegative  integer  and 


z$f(  r)  =  rL  z^(e,(p)  =  (-1)M  rL  P^icosd)  cosM(P  0  <M<L 

z~lm{  r)  =  rL  z2m(9,  </>)  =  (-1)M  rL  P*f  (cos  9)  sin  M<\>  1  <M  <L  (6) 

where  P/^  (cos  9)  denotes  a  Legendre  function  (see  0  Eq  8.751.1).  Notice  that  index  n  is  not 
the  usual  N  quantum  number,  but  it  is  related  to  it  by:  n  =  N  —  L.  Index  n  is  useful  to  simplify 
the  notation  in  many  equations. 

STO  can  be  normalized  by  multiplying  by  the  norm  factor: 


Nn 

lyLM 


2L  +  1  (L- \M\)\  {2Q2n+2L+l 


(7) 


_  2tt(1  +  (5m,o)  (L  +  \M\)\  (2n  +  2L)\ 

Eqs(|T|  to  (|4])  show  that  molecular  integrals  depend  on  products  of  two  functions  rather  than  on 
functions  themselves.  These  products: 


DlML'M'  —  Xlm(C:ta)  Xi'M'(C)rs)  (8) 

will  be  called  charge  distributions  (or  simply  distributions )  in  the  sequel,  a  name  coming  from 
the  role  played  by  these  quantities  in  the  integrals. 

Details  on  properties  of  STO  and  their  charge  distributions  are  given  in  reference  [0. 


3  Ellipsoidal  coordinates 


As  it  is  well  known,  ellipsoidal  coordinates  are  often  a  good  choice  to  deal  with  systems  having 
cylindrical  symmetry.  Since  this  is  the  case  of  all  the  two-center  integrals,  several  algorithms 
developed  for  these  integrals  are  formulated  using  this  type  of  coordinates.  Here  we  summarize 
their  most  relevant  features. 

Working  on  a  lined-up  system,  ellipsoidal  coordinates  £,  rj  y  0  associated  to  a  pair  of  centers,  A 
and  B,  lying  on  the  Z  axis  are  related  to  spherical  coordinates  by: 


ta  +  rB 

R 


V 


ta  ~  rB 

R 


(9) 
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and 


rA  =  (t  +  V)^ 


rB  =  (£-??)§ 


(f  +  v) 


(f  -  v) 


cos  9  a  = 


£v  + 1 

£  +  V 


COS0R  = 


£v  - 1 


It  must  be  noticed  also  that: 


R 


x 


=  xA  =  XB  =  —  a/(^2  -!)(!-  V2)  cos  < 


R 


y  =  VA  =  VB  =  j  y/ (£2  -  !)  (!  -  V2)  sin < 


R 

z  =  ZA  =  J  (€v  +  !) 

zB  =  ^  (tv  ~  !) 

In  this  system,  0  a  =  4>b  =  0,  and  the  volume  element  is  given  by: 


R3 


dr  =  —  (£2  -  rf)  drj  d£  dc^ 


and  the  definition  domain  is: 


(10) 


(11) 

(12) 

(13) 

(14) 

(15) 


0  <  0  <  2n 

—l<rj<l  (16) 

1  <  £  <  oo 

The  factor  1/ 1 r  —  r'|  can  be  expressed  in  terms  of  these  coordinates  as: 


r  —  r 


OO  l 


/  4  7  y(2  —  5m, o)  (2(  +  1) 

Z=0  m= 0 


(l  —  m)\ 
(. I  +  m)\ 


-i  2 


pr(t<)  q?(& 


x  Pzm(r/)  P™(t)')  [cos(m0)  cos(m0/)  +  sen(m0)  sen(m0')] 


(17) 


where  £<  =  min(£,  £'),  £>  =  max(00)  and  Pim(z),  Q)n(z)  are  the  corresponding  Legendre 
functions  -see  m  ,  seccion  8.7. 
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4  Translation  of  STO 


Molecular  integrals  are  specially  simple  if  all  the  functions  are  centered  at  the  same  point  of 
space.  This  fact  suggests  that  the  solution  of  integrals  involving  functions  centered  at  different 
points  could  be  simplified  by  referring  these  functions  to  a  common  center.  This  is  the  aim  of 
the  so-called  translation  methods,  which  are  actually  expansions  of  basis  functions  in  regular 
spherical  harmonics  times  radial  factors.  Placing  the  origin  of  coordinates  at  the  common  center, 
they  read: 


OO  l 

Xlm(C,ti)  =  ^m(r)  fL(r)  (18) 

1=0  m=—l 

and  the  several  translation  methods  basically  differ  in  the  particular  representation  chosen  for 
the  radial  factors. 

Barnett  and  Coulsonlfl2l  IT3ll  and  Lowdin|fl4l  IT5TI  proposed  more  than  fifty  years  ago  to  obtain 
the  radial  factors  from  the  addition  theorem  of  Bessel  functions  (Gegenbauer  cxpansionliTbl  ) 
which,  when  R/  lies  on  the  z  axis,  reads: 


where 


e~<ri 


f  i 


^2(21  +  1)  Pi  (cos  0) 
1=0 


(r  R7)V2 


(19) 


Mu((r,  C R)  =  I„(Cr<)  K((r>)  (20) 

r<  =  min (r,R),  r>  =  max(r,  R),  Pi(cos6)  are  the  Legendre  polynomials  (see  0  seq  8.91), 
0  is  the  angle  between  r  and  R/,  and  Iv{z),  Kv(z )  are  the  corresponding  Bessel  functions,  the 
latter  also  called  Macdonald  function  (see  0  sec  8.4-8. 5).  This  expression  is  easily  extended 
to  arbitrary  R/  bearing  in  mind  the  addition  theorem  of  the  Legendre  functions,  and  the  result 
is: 


e~Cri 


ri 


OO  l 

E  E  *r(Rj) 

1=0  m=—l 


(2 1  +  1)  (2  —  5m o)  ( l 
( l  +  \m\)\ 


m\)\Ml+1/2(Cr,CRi) 

(r  R/)m/2 


(21) 


This  equation  is  sufficient  to  study  the  simplest  integrals  needed  in  the  shift-operator  context 
to  be  described  in  the  next  section  but,  in  a  different  context,  it  must  be  extended  to  integrals 
involving  functions  with  higher  quantum  numbers.  An  efficient  way  to  carry  out  this  task  is  to 
combine  Steinbom’s  formulalITTll : 


d  CrR 

-  ^M„«r,CR)  =  [Af„_,(Cr,Cfi)  -  M„+l((r,(R)}  (22) 

with  the  recurrence  relationllT8l: 
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r2  +  R2  + 


(C  rRf 


2(i/ -  l)(i/  +  l). 


Mv(Qr\  (R)  =  rR 


v 


1  il4_i(Cr,  C-R)  +  — — -^i/+1(Cr;  CR) 


z/ 


CrR 


il4_2(Cr,CR)  +  Mu+2((rXR) 


(23) 


z/(z/  —  1)  z/(z/  +  1) 

for  increasing  the  n  quantum  number  and,  next,  to  use  the  recurrence  relations  of  the  Legendre 
functions  for  increasing  l  and  m.  The  compact  results  obtained  in  this  way  are  illustrated  with 
the  examples  of  appendix  1  of  refll20ll. 

The  radial  factors  of  eq(  2 1 )  and  their  generalizations  have  different  definitions  for  r  <  Rj  and 
r  >  Rj,  a  fact  that  can  be  unsuitable  for  some  applications.  The  alternative  are  the  one-range 
formulas,  in  particular  the  remarkably  simple  one: 


e-CrB 


TB 


2C  e“c  (R+r)  5^(2/  +  1)  Pi( cos  0)  (4C2  r  R)1 
1=0 


X 


oo 


E 


p\ 

( p  +  l+  1  )(p  +  21  +  1)! 


Lpl+1(2(  r)  L2J+1(2(R) 


(24) 


where  L2l+1(z )  is  a  Laguerre  polynomial  (see  0  seq  8.97).  This  expression,  which  was  first 
derived  by  usll2Tll22Tl  and  next  studied  by  Stcinborn[[23l,  is  valid  for  R/  lying  on  the  z  axis,  but 
can  be  generalized  to  any  orientation  of  R/  and  any  quantum  numbers  like  described  in  the  two- 
range  case.  There  exist  other  alternatives  obtained  by  projection  of  the  radial  functions [24,  25], 
but  they  are  more  complicated  than  eq(|24|)  and  its  generalization. 


5  The  shift-operator  technique 

A  very  powerful  tool  for  dealing  with  two-center  integrals  has  been  developed  in  our  group[[8lj9, 
IH,  and  it  is  based  on  the  fact  that  every  basis  function  can  be  obtained  by  applying  a  differential 
operator  which  acts  on  the  real  parameters  of  the  simplest  function.  For  STO,  it  gives: 

xL(C,r /)  =  nsB(j)^  (25) 

where 


nwn  =  *r(v,)(-D"  (-!!)'  <26) 

with  r/  =  r  —  R/,  and  2j"(Vj)  being  the  operator  obtained  by  replacing  the  corresponding 
cartesian  coordinates  by  the  derivatives  with  respect  to  them,  in  the  expression  of  the  regular 
harmonics  in  terms  of  these  coordinates. 
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Note  that  the  shift-operator  of  eqs(  25 )  and  (26)  relate  an  arbitrary  STO  with  the  simplest  one 
through  parameters  that  remain  after  integration  over  the  space  variables,  and  therefore  they 
also  relate  integrals  involving  arbitrary  STO  with  those  involving  the  simplest  STO.  As  a  con¬ 
sequence,  integrals  can  be  obtained  in  two  separated  steps:  first,  the  expression  of  the  integral 
involving  the  simplest  STO  is  derived,  and  next  the  shift  operator  is  applied  to  obtain  the  final 
integral.  This  approach  has  been  applied  to  several  types  of  integrals,  as  it  will  be  shown  later. 
The  master  formula  for  a  two-center  integral,  f^rL]J'r' ,  involving  arbitrary  quantum  numbers  can 
be  derived  from  the  basic  integral,  by: 


rn'L'M' 

JnLM 


^ I'M 


000 


d_ 

dC 


'{J)  /ooo 

d 

d~C 


-A77 


X 


rOOO 

J  000 


(27) 


where  the  basic  integral  f ^  =  f(R,(,C )  is  the  simplest  of  the  type  under  consideration, 
i.e.,  that  involving  two  0s  functions.  In  this  expression,  (  and  ('  are  the  exponents  of  the  0s 
functions,  and  R  =  X  i  +  Y  j  +  Z  k  =  (Xj  -  Xr)  i  +  (Yj  —  Yr)  j  +  {Zj-  Z »  k. 

Bearing  in  mind  Hobson’s  theoremlfTTfl  and  the  relation: 


d  d  d 

dX  =  dX]  =  _ dX j 

and  likewise  for  ^  y  ^ : 


(28) 


(V/)  z%'{Vj)  f(R,C,C)  = 


x 


L< 

(- i)L  E 


k= 0 


2~k 

~k\ 


n  o  \L+L'~k 

\RdRj 


f(R,CC) 


where  L<  =  min(L,  L'). 

Expanding  in  harmonics  the  product  of  regular  spherical  harmonics: 


(29) 


ZM 


:r)^"'(r)  =  EE' 


LML'M' 


'L-\-L'—cllm  “ L-\-L 


~L+U- 2,(R)  R 


21 


(30) 


1=0  m 


and  applying  V2  to  both  sides  of  the  equation,  the  first  term  of  29  is  attained.  Applying  once,  it 
leads  to: 


V2  4n(R)  R21  =  4 1  (p  +  l  +  1/2)  ^(R)  R2l~2 


(31) 
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and  iterating: 


V2fc  z™(R)  R 21 


and,  therefore: 


22k  l\  rQo  +  Z  +  3/2) 
(l-k)\  r(p  +  l  -  k  +  3/2) 


2™(R)  R2l~2k 


(32) 


(-1)L 

k= 0 


/  1  d\  L+L'~k 

\RdRj 


f(R,C,C)  (33) 


where: 


gpLML'M 


2^  A  /!  T(L  +  L'  -l  +  3/2)  R2f  2fc  lml’m'  m  /R'1 

k\  ^  (l  —  k)\  T(L  +  L1  —  l  —  k  +  3/2)  ^  L+L'-2im  l+l'-iA  ) 

l=k  '  m 

(34) 


aL+L'-2im  being  the  expansion  coefficients  of  the  product  of  regular  spherical  harmonics  in 


harmonics. 


with: 


rn'L'M'  _ 
inLM  ~ 


=  (-l) 


L< 

LE 

k=0 


'(R )  fZLn'L'(R, (,(') 


(35) 


fl 


■«**  -  i-sr  (-* 


1  d_ 
'C  5C 


I A 

C  5C7 


X 


1  d  \  L+L'~k 

urn) 


(36) 


The  presence  in  eq(  36 )  of  Bessel  and  derivative  operators,  which  do  not  commute  with  one 
another,  can  complicate  the  derivation  of  the  master  formula,  but  this  can  be  circumvented  by 
using  the  following  identities: 


or 


1  d_ 

'(  dC 


E(n/  2) 

:-i)”+I  c  E 


n\ 


£ 


n\  (2C) 


2=0 
2  i—n 


(n  -  2i)!  i!  (2C2)* 


(n  —  *)!  (2i  —  n)!  2*  <9£ 


1  <9 


i-\-L 


dy-i+L 

U  d~c) 


(37) 
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nL 


=  (i  +  P-  1)!  E 


n! 


(L  —  i  —  1)!  (n  —  p  +  i)\  (p  —  i)\  2i  i\ 


Alternatively: 


with 


(39) 


(40) 


rnL  / 

cp  =  (n 


+  2L  —  p  —  2)\  ]T 


n\ 


(41) 


(L  —  i  —  1)!  (n  —  p  +  i)\  {p  —  i)\  2L~1~l  i\ 

In  both  (39)  and  (|4T|)  the  summation  runs  over  the  positive  values  of  i  for  which  the  factorials 
make  sense. 

Obviously,  the  final  formula  and  the  algorithm  stability  will  depend  on  the  expression  taken  for 
the  basic  integral  f(R,  ('). 


6  Two-center  one-electron  integrals 

These  integrals  can  be  classified  in  three  groups:  two-center  overlap: 


,  Xlm  I  Xl'm'  )  —  dr  X'lat(C?  Xl'm'(C  >  rs) 


kinetic  energy: 


(42) 


(  Xlm  I  —  2^2  I  X^vm'  )  —  Xlm(C,  rA)X72XL'M'(C,  rs)  (43) 


and  nuclear  attraction  integrals: 


(  Xlm  I  I  Xl'M'  )  =  ~Qa  f  dr  Xlm((i  fa)  Xl'M'((' j  rs)  (44) 

rA  J  rA 


The  kinetic  energy  and  nuclear  attraction  integrals  of  eqs([43j)  and  (|44|)  can  be  easily  expressed 
in  terms  of  overlap  integrals.  In  fact,  since: 
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V^M-K',  r)  =  C'2  xiwi C',  r)  -  2C'  (£'  +  »')  x£v‘.(C',  r)  +  (rf  - 1)  (2 L'  +  n')  x^(C', r) 

(45) 

it  follows: 


(  Xlm 


-  Y  (  Xlm  I  xf-M' )  +  C'  (£'  +  n')  ( xnLM  I  X%-M\  ) 
in'  -  1)  (L  +  n'/2)  (  X"lm  I  x£„2  )  (46) 


Moreover,  it  is  also  obvious  that: 


(  Xlm  I  I  Xz/m'  )  ~  Qa  (  Xlm  I  Xl/m'  )  (47) 

rA 

Clearly,  the  analysis  of  the  computation  of  the  two-center  one-electron  integrals  can  be  reduced 
to  the  overlap  integrals.  These  are,  in  principle,  complicated  functions  depending  on  five  real 
variables  ((,  XB  —  Xa,  Yb  —  Ya ,  ZB  —  Za )  and  six  integer  quantum  numbers  (n,  L,  M, 
n',  L',  M').  However,  using  normalized  functions  and  taking  a  lined-up  axis  system,  they  are 
reduced  to  functions  of  two  real  variables  ((  RAB ,  ('  Rab )  and  five  indices:  (n,  L,  n' ,  L',  M ) 
which  can  be  further  reduced  to  two-indices  functions  depending  on  two-variables  by  means  of 
the  recurrence  relations  of  the  spherical  harmonics. 

Simple  and  appealing  algorithms  for  their  calculation  can  be  designed  by  reversing  this  process, 
i.e.  by  computing  first  the  two-indices  functions  and  applying  next  the  recurrence  relations  to 
increase  the  quantum  numbers.  Notice  however  that  there  are  different  possible  choices  for  this 
pair  of  indices,  and  every  choice  determines  the  set  of  recurrence  relations  to  be  applied.  In  prin¬ 
ciple,  the  best  choice  should  lead  to  (i)  accurate  and  nonexpensive  two-indices  functions,  and 
(ii)  a  stable  recursion  scheme.  Examples  of  two-indices  functions  fulfilling  the  first  requirement 
are: 


and 


(48) 


( Xoo  Xm }  =  J  dr  r"-1  t%-1  (49) 

The  first  one,  eq(|48]),  was  the  choice  in  the  first  version  of  SMILES,  in  spite  of  the  fact  that 
the  associated  recursions  are  not  stable,  because  that  version  was  intended  for  relatively  low 
quantum  numbers,  in  which  case  the  recurrence  relations  had  to  be  applied  a  rather  reduced 
number  of  times. 

The  second  choice,  eq(|49]),  is  that  implemented  in  the  current  version  of  SMILES.  It  has  two  im¬ 
portant  features:  the  starting  two-indices  integrals  can  be  computed  with  full  accuracy,  and  the 
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associated  recursions,  though  not  being  fully  stable,  downgrades  more  slowly  than  the  previous 
one,  so  that  moderate  quantum  numbers  can  be  reached  with  a  reasonable  accuracy. 
Nonetheless,  as  it  will  be  shown  in  the  next  section,  it  is  still  not  sufficiently  good  for  very 
high  quantum  numbers  and,  therefore,  we  have  developed  a  new  algorithm  which  improves  the 
performance  of  the  previous  ones. 

6.1  Algorithm  1:  the  old  scheme  (( Xoo  I  Xoo  ))  scheme 

The  procedure  starts  with  the  two-indices  integrals  ( Xoo  |  Xoo  )  'n  the  lined-up  axis  system,  and 
next  uses  recurrence  relations  for  increasing  the  remaining  quantum  numbers.  The  starting 
integrals  are: 


(50) 


Different  algorithms [U  lead  to  different  expressions  for  the  pending  integral,  and  we  must 
choose  the  most  convenient  one  for  taking  the  derivatives  of  eq(|50|)  and  preserving  full  ac¬ 
curacy.  This  is: 


( Xoo  I  Xoo)  =  /  dr 


e~Cr  e-CrB 

r  rB 


4  7T 

c  +  c 


du  e~Cru  e-f'  *(!-«) 


(51) 


B 


where  R  is  the  distance  AB:  R  =  |Rg  —  R^|  =  |Rg|  =  | Z 
Replacing  ([5T])  in  (50),  the  derivatives  can  be  easily  taken  with  the  aid  of  Leibnitz  rule,  leading 
to: 


[  Xoo  I  Xoo  )  = 


p=0  p'= o 


n 


P )  \P 


n'\  (n  +  n'  —  p  —  p')\  Rp+p' 


(C  +  CO 


n+n'—  p— p'+l 


-1 


—  P~iRu  P-CR(i-u) 


x  /  du  up  (1  —  u)p  e  e 

Jo 


(52) 


or  equivalently: 


1 

(n  +  n1  +  1  —  k)\  [(£  +  (')  R]k 
X  (h  -p)\ p\  n '  +  1  +  P  -  fc;  n  +  n'  +  2  -  fc;  (C  -  C')/R)  C  >  C  (53) 

p=max(0, /c—n) 

where  0(a,  /3,  z)  is  the  confluent  hypergeometric  series  (eq  9.210.1  of  reflTTlI). 


,  „  .  /  4  7T  n\  n'\  Rn+n'  e  R 

,  Xoo  I  Xoo  )  = - -  2^ 

k= 0 


C  +  C 
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Eqs(52)  and  (|53|)  keep  the  whole  accuracy  if  (  >  ('  but  fail  for  (  <  In  this  case,  accuracy  is 
recovered  by  means  of  Kummer’s  transformation: 


0M,  (C  -  CO  R)  =  e-«~VR  </>(/?  -  a,/?,  (C'  -  C)  R)  (54) 

which  simply  exchanges  the  order  in  the  pairs  (£,  (')  and  (n.  n'). 

The  set  of  recurrence  relations  needed  for  increasing  the  remaining  quantum  numbers  can  be 
easily  derived  from  the  definition  (M  >  0): 

z%( r)  =  {M  (2  r  sin0)M  cos  M<t>  (55) 

v  7r 

and 


z 


M 
L+ 1 


1 

L-M  +  l 


(2 L  +  1)  r  cos 6  Zjf (r)  -  (l  +  M)  r2  z?_ i(r) 


with  the  convention:  z}f(r)  =  0  for  L  <  M. 

Combination  of  eq([55|)  with  the  cosine  theorem  gives  the  first  recurrence  relation: 


(56) 


(  Xm+im+i  I  Xm+im+i  )  —  ~  [^(  Xmm  I  Xmm  )  +  2i?2(  Xm~m  I  Xmm  ) 


+  2 R2(Xmm  I  Am  )  (  Xmm  I  Xmm  )  (  Xmm  I  Xmm  )  7?4(  Xmm  I  Xmm  )]  (57) 

and  eq(f56])  leads  to: 


1  xn'  )  -  2L+l 

(  Xlm 

Xl'm  )  +  R2  (  Xlm 

Xl'm  )  (  Xlm 

X 

^  + 

^  to 

\Xum)  L-\M\  +  1 

2  R 

(L  +  \M\)  (x1-im  I  Xlm) 


(58) 


1  \n'  )  -  2L  +  1 

r  (xiti 

Xl'm  )  R2  (  Xlm 

Xl’m  )  (  Xlm 

n’+2  \ 

XVM  / 

IXv+iM)  L-\M\  +  1 

2  R 

(L+|M|)(x2mIxS21m) 


(59) 


The  recursion  scheme  proceeds  as  it  follows: 


M  = 

0  : 

( Xoo 

X  oo ) 

|58l> 

© 

M  = 

1  : 

<xj, 

Xnn) 

(58) 

© 

M  = 

2  : 

Xlo 

Xoo  ) 

([59l> 

<  Xlo 

Xl’o  ) 

Xli 

\xnu) 

j59} 

(xl, 

1  xin ) 

(60) 
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6.2  Algorithm  2:  the  ellipsoidal  coordinates  scheme 

A  new  algorithm  developed  in  this  project  is  based  on  the  formulation  of  the  problem  in  ellip¬ 
soidal  coordinates.  This  coordinates  system  implies  working  in  a  lined-up  system  as  discussed 
in  the  previous  section  in  which  integrals  involving  funtions  with  different  values  of  the  M 
quantum  number  are  null  by  symmetry.  Therefore,  we  will  consider  only  the  case: 


,  Xlm 


Xl'm  )  =  drz?{rA)z%(rB) 


„n-l  n'-l  p-(rA  -C,'  rB 
A  rB  e  e 


(61) 


To  express  these  integrals  in  ellipsoidal  coordinates,  we  will  separate  the  representation  of  the 
spherical  harmonics  from  that  of  the  radial  factors.  For  these  latter,  it  reads: 


(p  \  n+n'+l 

-  )  (f  +  v)n  (C  -  n)n'  e~K  e~vr>  (62) 

where  ft  =  ( C  +  CO  -R/2  and  u  =  {Q  —  Qr)  Rj 2.  Expanding  the  binomials  and  grouping  terms, 
it  results: 


dr  r™”1  rnB-x  e~<rA  e~c'rs  =  d£  dr]  d(f)  fj) 


n+n'+l 


e~K  e-"v 


n+n' 

in+n'~P  rf  a™'  (63) 

p= o 


with 


min  (n,p) 


C'  =  (-i)? 


n\  n'\  (— 1)* 


^  (n  —  i)\  (n  —  p  +  i)\  (p  —  i)\  i\ 

i=max(0 ,p-ri)  K  y  ^  ’ 


The  coefficients  ar"1  can  be  recursively  obtained  by  means  of: 


and 


n+ln'  _  nn'  ,  nn' 

ap  ~  ap  '  up-l 


nn! 1  _  nn'  _  nn' 

ap  —  up  up_1 


Eq(|63])  leads  to: 


(j~ j  \  1 

Y,  An+„'-M  B„A) 


p= 0 


where 


(64) 


(65) 


(66) 


(67) 


Aj((3)  =  J  d££J  e 


3 


(68) 
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and 


B,j(u)  —  /  dr/ rf  e 


—  V  T) 


'-1 


Let  us  consider  now  the  expansion  of  the  spherical  harmonics.  Since,  (  Xlm  I  Xl'm  )  = 
we  will  consider  only  the  case  M  >  0,  in  which: 


(69) 

(  Xl-M  I  Xl1  -M  ) 


(fa)  M  =  2  (!  +  cos  2M0)  -  1)"  (1  -  tf) 


2 \M 


X 


X 


l-m\ 

1  VA  '  (L-l/2-i)!  (-1)' 


7T 


E 

i=0 


^v+l)L-M-2i  (e  +  r?) 


(L  —  M  —  2i)\  i\  22i 


\2  i 


(  L'  —  M  \ 

1_  '  (1/  —  1/2  —  i')\  (-1) 

i'= o  v  ' 


Expanding  in  powers  of  £  and  //  and  grouping  terms,  it  follows: 

L+L’  L+L'  2k 


c“(r.- 1)  z“(rs)  =  i  (1  +  cos  2 M</>)  (  ^ 


j2J2ek~w  no 


fc=0  j=0 

Multiplying  ([63])  by  (|7T|),  the  general  integral  can  be  expressed  as: 


(  xIm  I  Xu+1M  )  =  tt  (1  +  5M0)  fj) 


n+n'+L+L'+l  L+L'  2k  n+ri 

EEE  -/4n+n/+2fc— j— p  (P) 

k= 0  j= 0  p=0 


v  R  „LML'M 

X  Bj+p(V)  Cl  Olkj 


(72) 


or,  alternatively: 


(XlmIXl'm)  —  /(  (1  +  ^mo)  (  2 


/?, 


7T-— |— 77-/  — [— X/— (— — |—  1  L~\~Lf  TL-^Tl' ~\~2k 


X 


E  E  ■A-n+n' +2k— A  (/7)  SA(i/) 

fc=0  A=0 

min(A,2fc) 

E  nn'  LML'M 

°A -j  akj 

j=m&x(0,\—n—n') 

/  \  n-\-nf +1  L-\-L r  n-\-n' -\-e2k 

7T  (1  +  SmO )  (  ^ 


E  E  A-n-\-n' +2k— A  (/3)  Sa(iz) 


fc=0  A=0 


X 


min(A,A— ?i— n') 

,.nn'  „  LML’M 

p=max(0,A —n—n') 


ap  ak\-p 


(73) 
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The  coefficients  a^fL'M  can  be  obtained  by  recursion.  We  start  from  eq(  7 1 )  for  the  particular 
case  L  =  LI  =  AT : 


M/  \  Mi 


and  define: 


1 

2 

1 

2 


(  r\2M 

i 

? 

1 

to 

_ 1 

V  2  / 

fn 

(  R  \2M 

I 

to 

V  2  / 

2m  (f  -  1)M  (1  -  if) 

2  2  M  2k 


2 \M 


2" 


2^  2 j  ^2i 


fc=0  *=0 


X 


(-l)fc+M  (AT!)5 


(AT  +  fc  +  i)!  (A7  -i)!  (Jfc-i)!  i! 


(74) 


„ MMMM  n 
ak2i+l  =  0 

(_1  )k+M  22M  [(AT  -  1/2)!]2  (AT!)2 
7T  (AT  —  A;  +  i)!  (AT  —  i)\  ( k  —  i)\  i\ 

with  0  <  i  <  k. 

Next,  the  recurrence  relation  for  increasing  L: 


„  MMMM 
ak2i 


(75) 

(76) 


*L+ i(f) 


leads  to: 


1 


L  -  AT  +  1 


L+1ML'M 

akj 


where  z^_f  r)  =  0. 
For  increasing  L 


zu+i  (r) 


leads  to: 


L'-M  +  l 


(2L+1)  |  (l+£„)  zf(r)-(L+M)  (|)  K2+2£„+„2)  z£i(r) 

(77) 


277  +  1  /  LML'M  |  LML'M 

L-M  + 1  ^ 

L  +  AT  /  L—1ML'M  ,  o  „L-IML'M  ,  L—IML'M 

L-M  +  l  \  k~^  +2ak-ij-i  +ak- ii-2 


(78) 


— (2Z/+1)  |  (1-^)  z"(r)-(L'+M)  (+)"  K2-2eW) 

(79) 


LML'+IM 

akj 


2L'  +  1 


a 


LML'M  „  LML’M 


kj 


L'-M+l 
L  +  M  /  lml'-im 


—  a 


L'-M  +  l 


fc-i j-i 


r,  ^LML'-IM  t  „LML'-1M 
~  2  afc_i  j_i  +afc_li_2 


(80) 
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6.3  Tests  on  the  accuracy  of  algorithms  1  and  2 

The  first  test  illustrates  the  accuracy  of  the  two  above  described  algorithms  in  case  of  moderate 
quantum  numbers.  Fig  |T|  gives  the  order  of  magnitude  of  the  integrals  involving  normalized 
STO  for  a  range  of  exponents  and  with  quantum  numbers:  NA  =  NB  =  10,  LA  =  LB  =  5, 
and  —  5  <  M  <  5.  Notice  that,  in  the  notation  used  throughout  this  report:  n  =  N  —  L,  i.e., 
the  integrals  correspond  to  nA  =  nB  =  5  in  this  notation.  The  scaled  exponents  (£A  =  (A  R, 
C B  =  (B  R )  range  from  5  •  10  :!  to  5  ■  102  in  an  almost  exponential  way. 


Integrals  n  =n  =5;  I  =1=5;  R  =  1 

a  A  B  ’  A  B  ’ 

Overlap  integrals  with  normalized  STO:  nA=n=5, 1=1=5,  R  =  1 


Decimal  logarithm  of  the  absolute  value  of  the  integrals 
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Figure  1:  Overlap  integrals  with  NA  =  NB  =  10,  LA  =  LB  =  5 

Region  of  exponents  corresponding  to  negligible  integrals  (absolute  value  lower  than  10” 15)  is 
coloured  in  this  and  the  following  figures.  The  behavior  of  the  algorithms  in  this  region  is  rather 
irrelevant,  but  we  have  kept  the  values  attained  therein  to  fully  illustrate  the  analysis  carried  out. 
It  must  be  recalled  that  the  information  on  the  number  of  absolute  decimal  figures  in  this  region 
can  be  completely  useless:  for  very  small  integrals,  there  can  be  a  high  number  of  accurate 
decimal  figures  (zeroes)  even  in  cases  in  which  the  order  of  magnitude  in  the  computed  integral 
is  wrong. 

The  figure  shows  that  these  integrals  have  meaningful  values  even  for  relatively  large  scaled 
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exponents  (in  the  range  of  30). 

The  number  of  correct  significant  figures  (i.e.  the  relative  error)  attained  with  both  algorithms 
working  in  double  precision  are  reported  in  figure  |2j  In  this  case  as  well  as  in  the  remaining 
tables,  the  reference  values  were  computed  with  the  algorithm  based  in  ellipsoidal  coordinates 
using  multiprecision  (65  significant  digits  in  computation).  A  previous  analysis  on  the  consis¬ 
tence  of  the  results  was  made  using  different  working  precision  values  (lower  and  higher  that 
those  finally  chosen)  in  a  number  of  selected  cases.  This  analysis  showed  that  the  number  of 
digits  finally  taken  in  the  calculations  was  sufficient  to  yield  more  than  20  correct  significant 
figures  in  the  results  for  all  the  selected  range  of  exponents.  Figure |3]gives  the  number  of  correct 
decimal  places  (absolute  error). 
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Integrals  n=n=5;  I  =1=5;  R  =  1 

a  A  B  ’  A  B  ’ 

Algorithm  based  on  recursion  over  L,L'  and  M  starting  from  integrals  <N00|N'00> 


Number  of  accurate  significant  figures  of  integrals  computed  in  double  precision 
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Integrals  n  =n  =5;  I  =1  =5;  R  =  1 

a  A  B  ’  A  B  ’ 

Algorithm  based  on  ellipsoidal  coordinates 


Number  of  accurate  significant  figures  of  integrals  computed  in  double  precision 
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Figure  2:  Accurate  significant  figures  in  double  precision  for  overlap  integrals  with  Na  =  Nb  =  10,  La  =  Lb 
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Integrals  n=n=5;  I  =1=5;  R  =  1 

a  A  B  ’  A  B  ’ 

Algorithm  based  on  recursion  over  L,L'  and  M  starting  from  integrals  <N00|N'00> 


Number  of  accurate  decimal  figures  of  integrals  computed  in  double  precision 
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Integrals  n  =n  =5;  I  =1  =5;  R  =  1 

a  A  B  ’  A  B  ’ 

Algorithm  based  on  ellipsoidal  coordinates 


Number  of  accurate  decimal  figures  of  integrals  computed  in  double  precision 
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Figure  3:  Accurate  decimal  figures  in  double  precision  for  overlap  integrals  with  N_ 4  =  Nb  =  10,  La  =  Lb  =  5 
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Integrals  n=n=5;  I  =1=5;  R  =  1 

a  A  B  ’  A  B  ’ 

Algorithm  based  on  recursion  over  L,L'  and  M  starting  from  integrals  <N00|N'00> 


Number  of  accurate  significant  figures  of  integrals  computed  in  quadruple  precision 
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Figure  4:  Accurate  significant  figures  in  quadruple  precision  for  overlap  integrals  with  Na  =  Nb  =  10,  La 

Lb  =  5 
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Figure  5:  Accurate  decimal  figures  in  quadruple  precision  for  overlap  integrals  with  Na  =  Nb  =  10,  La  =  Lb 
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Figure  6:  Overlap  integrals  with  Na  =  Ng  =  20,  L  A  =  LB  =  10 
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Figure  7:  Accurate  significant  figures  in  double  precision  for  overlap  integrals  with  Na  =  Nb  =  20,  La  =  Lb 
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Figure  8:  Accurate  decimal  figures  in  double  precision  for  overlap  integrals  with  Na  =  Nb  =  20,  La  L />  =  1 0 
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Figure  9:  Accurate  significant  figures  in  quadruple  precision  for  overlap  integrals  with  Na  =  Nb  =  20,  La 
Lb  =  10 
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Figure  10:  Accurate  decimal  figures  in  quadruple  precision  for  overlap  integrals  with  Na  =  Nb  -  20,  La 
Lb  =  10 

From  these  figures,  it  is  clear  that: 
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•  Algorithms  1  and  2  are  complementary  in  the  sense  that  algorithm  1  works  better  for  large 
scaled  exponents  whilst  algorithm  2  works  better  for  smaller  exponents. 

•  Though  the  accuracy  of  algorithm  2  is  roughly  higher,  neither  algorithm  1  or  2  in  double 
precision  give  sufficient  accuracy  in  the  whole  range  of  exponents. 

Figures  [4]  and  [5] enable  the  same  comparison  in  case  of  integrals  computed  in  quadruple  preci¬ 
sion.  It  is  clear  now  that 

•  Algorithm  1  can  be  confidently  used  except  for  very  small  scaled  exponents  (lower  than 
0.15)  whilst  algorithm  2  gives  satisfactory  results  in  the  whole  range  of  exponents  studied. 

•  Algorithm  2  downgrades  when  one  of  the  exponents  is  very  high  (~  500)  and  the  other, 
moderate  (~  8). 

In  the  second  test,  the  results  of  both  algorithms  are  compared  in  a  case  of  higher  quantum 
numbers:  Na  =  NB  =  20,  La  =  LB  =  10,  —10  <  M  <  10  (n^  =  nB  =  10  in  our  notation). 
Figure  [6]  shows  the  order  of  magnitude  of  these  integrals  as  a  function  of  the  scaled  exponents 
in  the  same  range  as  before.  Now  integrals  have  meaningful  values  for  exponents  as  large  as 
60,  notably  higher  than  in  test  1 . 

Figs  [7]  and  [8]  fully  confirm  the  previous  conclusions  reached  in  test  1:  algorithms  are  comple¬ 
mentary  and  double  precision  is  not  sufficient.  Furthermore,  comparing  figs  [2] and  [3]  with  [7] and 
[8]  it  becomes  evident  the  quick  downgrading  of  the  results  of  algorithm  1  in  double  precision  as 
the  La  and  LB  quantum  number  increase.  Figures  [9] and [T0| show  that  quadruple  precision  is  not 
sufficient  to  reach  an  acceptable  final  precision  in  these  integrals.  On  the  contrary,  algorithm  2 
still  works  well  in  almost  the  whole  range  of  scaled  exponents;  it  only  fails  when  an  exponent  is 
very  high  (~  500)  and  the  other  is  moderate  (~  8).  Fortunately,  these  integrals  are  small  (10~10 
to  10~14)  and  just  in  these  cases  algorithm  1  gives  sufficent  accuracy. 

The  algorithms  used  for  computing  the  Aj(/3)  and  Bj(v)  integrals  are  fully  stable  (see  appendix 
A),  and  so  are  the  computations  of  the  a™n'  and  ay  L'M  constants.  The  accuracy  losses  come 
from  the  cancellations  between  the  positive  and  negative  terms  appearing  in  the  summation  of 
eq(|73]). 


6.4  Shift  operators  for  two-center  one-electron  integrals:  an  alternative  algorithm 


The  algorithm  described  in  the  previous  section  6.2  has  proved  to  be  sufficiently  accurate  for 
one-electron  two-center  integrals  in  a  wide  range  of  exponent  values  for  fairly  high  quantum 
numbers  and  it  has  been  adopted  for  computing  this  type  of  integrals.  Nevertheless,  as  some 
algorithms  for  Coulomb  integrals  that  will  be  described  below  require  overlap  integrals  with 
rather  high  quantum  numbers,  a  third  algorithm  has  been  developed  to  compute  them  in  these 
extreme  cases. 


31 


In  this  alternative  approach,  the  master  formula  for  the  overlap  integral  is  written  as: 


,  Xlm  I  Xl 


,'M'  )  ~  (“1)L  X/ 


yLML'M'  r^  \  gnLn'L' 


(R,C  C) 


k= 0 


where 


(81) 


^nL  n'  L 


'(R,CC)  = 


X 


I  JL\ 

RdRj 


f_d_\n'  fl  d_\L 

V  dCJ  VC  d() 

L+L'-k 

f(R,C,C) 


(82) 


Derivation  of  explicit  expressions  for  the  S\'.  L  n'L'  coefficients  is  simple  if  the  following  expres¬ 
sion  is  taken  for  the  basic  integral: 


( Os  |  Os' )  =  \Z8tt  R  /  duk_i/2((uR)  (83) 

Jo 

The  action  of  the  three  Bessel  operators  on  the  modified  Macdonald  functions  is  straightforward 
and  gives: 


S^L-'L'(R,CO  =  (~l)L+L'+t  Vstr  R2k+1  (-jf) 


x  /  du  uL  (1  -  u)L'  Lfc— 1/2(C«  R) 
Jo 


(84) 


There  are  now  two  possibilities.  If  the  derivative  operators  are  expressed  in  terms  of  Bessel 


operators  with  the  aid  of  eq(37 ),  after  some  algebra,  it  comes: 


n  n 


s^l\r,u')  =  (  E  E  ctm 


(2C)“  PC')”' 


£[  VeiI] 


x  I  du  uL+l  (1  -  u)L'+l'  /c_fc_i_i/_i/2(CUJR) 


L'+i' 


(85) 


where  0  <  k  <  min(L,  L')  and: 


c?m  = 


(_l)n+i  n|  (2(2  R2 y 


(86) 


(2 i  —  n)\  ( n  —  i)\ 

The  pending  integral  can  be  computed  by  numerical  methods.  Alternatively,  the  problem  can 
be  reduced  to  compute  integrals  with  the  general  form: 
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1/2(W’W')  =  I  dliv?  ( 1  -  u)P'  k-m- 1/2  (^\J  w2u  +  m/2(1  -  u)^j  (87) 
which  have  been  discussed  elsewhere  [fT8l  IT91. 

In  the  second  possibility,  the  integral  of  cq(f84])  can  be  solved  and,  next,  the  derivative  operators 
are  applied.  In  order  to  do  this,  we  carry  ut  the  expansion  of  the  Macdonald  function  in  Gegen- 
bauer  polynomials  and  integate  term  by  term.  After  some  cumbersome  algebra,  the  final  result 
can  be  written  as: 


gnLn'  V 


d  x  n 


d  y  (-i)^  L+^2k  LL, 

/  ,  ^k+l/2\J)  Mk+l/2+j 


d(' 


(V2  _  ^k+1/2 


3=0 


where 


^k+1/2 (j)  =  (-!)J 


x 


(k  +  1/2  +  j)  (2 k  +  j)\ 

2fc— 1/2  j\ 

(-l)i  (1/  —  k)\  (L  —  k  +  i)\  ( L  +  i)\ 
y  (1/  —  k  —  i)\  (L  —  k  +  i  —  j)\  (L  +  k  +  1  +  j  +  i)\  i\ 


L'-k 


£tT7 


and 


(88) 


(89) 


Mi 


k+l/2+j  =  h+l/2+j 


R{  C  -  CO 


K , 


2  i  --fc+i/2+j  i  2 

Defining  a  =  £  +  £'  and  a'  =  (  —  and  taking  into  account  that: 


R  (C  +  CO 


(90) 


d_\n 

^c;  V  aco 


n+n' 

:-i)n+”'  EC 

p=0 


d_ 

da 


n-\-n'—p 


_d_ 

da' 


with 


rnn'  _ 
LP 


min  (n,p) 

(-i)”  E 


(—1)*  n\  n'\ 


i=max(0, 
p—nf ) 


(n  —  i)\  ( n '  —  p  +  i)\  (p  —  i)\  i\ 


a  compact  expression  is  attained: 


(91) 


(92) 


S%Ln'L'  =  (_iy+n'+L+L'+k  ^  22k+l  /2(j) 


L+L'-2k 

E 

3=0 


n-\-n' 

x  E  <  F£V7/(n, K)  Gg1/2J(/, R) 

p= 0 


(93) 
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where: 


AS/2>.fl)  = 

/  0  \ 
\0a  / 

m  Afe+V2+>A/2) 

aAi+l/2 

(94) 

4+,/2+j(a'fl/2) 

(95) 

V  dcd 

)  afk+1'2 

Functions  F^^2j(a,  A )  and  G^_\^2  -{ot\  R)  have  different  expressions  depending  on  the  for¬ 
mula  taken  for  the  products  of  powers  by  Bessel  functions.  Thus,  taking  derivatives  directly 


with  the  aid  of  Leibnitz  rule  gives: 


(-1)’ 


and 


G 


(m) 

fc+1/2  ,j 


(k  -  1/2)!  am+k+1/2  \  i 


y:  |  m  j  (m  +  k  —  1/2  —  i)! 


x 


E 

r=0 


A' 


fc+l/2+j— i+2r 


(aR/2) 


(a',R)  = 


-l)r 


x 


(A: -1/2)!  cdm+fc+1/2  ^  V* 
E  0  -^fc+l/2+j— i+2r  (a'R/2) 


y]  (ra  +  /c  —  1/2  —  i)\ 


r=0 


a  R 


(96) 


a'  A 


(97) 


If,  as  an  aternative,  we  write: 


and 

„(m|  ,  ,  r,=  (R\MI2  (8_\m  (*_R V  WiW) 

11  ‘/2J  -  l  2 )  La7  V  2  )  (a'J?/2)‘+'/2+j 

Leibnitz  rule  gives: 


(98) 


(99) 


pM 
fe+l/2  j 


(cc,  A) 


\  fc+l/2+j 

2  J 


i  d  rv?-m 


m!  j !  cr 


E 

p=max(0, 

m-j) 


(-l)p  (aR/2)2P 
(m  —  p) !  (j  —  m  +  p) ! 


x 


\  ^  _ l-1) _ 

(p  —  2i)\  i\  (a2  A2/2)* 


k— k— 1/2— j—  p+i  (aR/2) 


(100) 
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and 


Kjrk+l/2,j 


(a',R) 


(R\k+1/2+j  ...  ,j-m  {a’ R/2fP 

V  2  /  m-J-a  S  (m-p)!  (j -m  +  p)! 

p=max(0, 

m-j) 


X 


£(§) 

1  Ik+l^+j+p-iioi'  R/2) 

^  (p-2i)!  i!  (a'2R2/2y  (a'  R/2)k+1/2+i+p~i 


(101) 


Notice  that,  whereas  eqs(|96])  and  (|100[)  are  free  of  cancellation  issues,  eqs(f97|)  and  (|99|)  may  be 
not.  Therefore,  the  former  are  preferable. 


7  Coulomb  integrals 

The  Coulomb  integrals  between  four  Slater  functions  are  defined  as: 


[xZma  Xl\m'a  I  XnLBMB  XvBM'B\  =  j  dr  I  dr>  xZmACa^a)  X^m'JCa^a) 

x  xTh  mb  (C b,  r's)  xLf  ML  (Cb,  r'B)  1_f,  (102) 


r  —  r 


Since  one-center  STO  distributions  can  be  expressed  in  terms  of  STO: 


xlm((a,y)  Xu  M'(c^r) — yy  yy 


LM  L  M  m  /  . 

aL+L'--2lm  ZL+L'-2l\r>  ' 


n+n'+2l-2  g-(C A+C,'A)r 


(103) 


l  m 


the  integrals  of  eq  ( 102 )  are  linear  combinations  of  more  simple  integrals  [Xlm  I  Xvm'  ]■ 


„  n  B 

Xlb  Me 


xd 


EEEE 


LaMa  L'aM'a  LbMb  L'bM'b  r  n  I  n'  1 
^  LA-\-L'A—2l  m  ^  LB+L'B—2l' m!  I  Xlm  I  A  I'm'  J 

(104) 


where  0  <  l  <  [(LA  +  L'A)/2],  0  <  V  <  [(Ls  +  L'B)/2],  each  the  sums  over  m  and  m!  contains 
two  terms  at  most,  and: 


[x?m  I  Xvm>  \  =  [Xim( C,  *a)  \  Xi'm'iC',  r'B )  ]  =  J  dr  J  dr'  x"m(C,  Xl>m'(C,  t'b) 
=  [dr  [  dr '  z™(rA)  rnAA~l  e~^A  zf( r'B)  r'Z”1  1 


|r  —  r'| 
(105) 


where  we  have  used  the  decomposition  of  the  product  of  regular  harmonics  in  harmonics,  and 
taken  (  =  (A  +  ('A  and  C'  =  (b  +  Ci b  ■ 

Eq|105|can  be  written  in  the  compact  form: 
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Xlm  I  Xl'vn' 


=  S„ 


dr  VJ”  (r)  xlC'M 


where 


dr  V,n(r )  |^m(r)  z?( rB)  r2k 


B 


3-C'  1-B 


rB 


(106) 


vr(r)  = 


An 


{n  +  2l  +  l)\  j{2l  +  l,(r) 


(2/  +  l)Cn+1l  (2/)! 

n—  1  r 

(Crf  £ 


—  e 


i= 0 


(C  r)2l+1 

(n  +  21  +  1)!  n\ 

(i  +  2l  +  1)!  ~~  ¥ 


(107) 


For  odd  n',e  ^'rB  /tb  must  be  replaced  by  e  ^  rs . 

A  first  general  expression  can  be  by  direct  substitution  of  the  potential  of  eq(]  107|)  into  the 


definition  of  eq(  106)  gives: 


Xlm  I  Xl’M'  ]  ~  Smm' 


An  f  (n  +  2L  +  1)! 


(2  L  +  1)  (  £n+2L+2 


dr 


Zl(* )  XnL'M’{rB) 
J.2L+1 


H  /**?( T)r 

3=0  ' 


j-2L-l  -C  r  n' 


e  Crx2'M'M 


n—  1 

E 


(n  +  2L  +  1)!  n! 

Cn+1  ^  L  (1  +  2L  +  1)!  _ 


0  /  drzf(r)rle  Cr  x2'M'M 


(108) 


The  last  term  in  the  r.h.s.  is  a  sum  of  ordinary  overlap  integrals  which  have  been  previously 
discussed.  The  first  term  can  be  regarded  as  an  overlap  integral  with  (  =  0  and  negative  index 
n.  The  second  term  is  a  sum  of  overlap  integrals  with  negative  n  index  that  can  be  attained 
either  by  translation  methods  or  by  recurrence  relations. 

This  integral  can  be  written  also  as: 


(n  +  2 L  +  1)!  (n'  +  2 L'  +  1)!  .  0  .  0 

(2 L  +  1)!  Cn  (2 L'  +  1)!  Cn'  XlM  Xl'm' 


An  (n  +  2L  +  1)! 

(2L'  +  1)  C/n'+1  (2L  +  1)!  Cn 


n! — 1 


E 


(n'  +  2L'  +  l)! 
(i  +  2  L'  +  1)! 


C  (  Xlm  I  Xz/m' 


) 


47T 

(2L  +  1)  C+1 


n—  1 


E 


(n  +  2L  +  l)! 
(i  +  2T  +  1)! 


(109) 
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A  second  general  formula  for  the  Coulomb  integral  of  eq(|106[)  is: 


r  ,  n  |  n'  1  r  a  _  \  ^  \  ^  r>LMkV  M  T}L+L'+2k—l—2p—l/2 

[  Xlm  I  XL’ M •  \  =  °MM’  4  7T  2^2^Bip  R 

l  p 

/»oo 

x  /  dr  V£(r)  r2p+l+'s/2  Ml+l/2  (110) 

Jo 

in  which  a  one-dimension  integral  appears,  that  can  be  solved  by  numerical  or  analytical  methods lfl~8l. 
As  an  alternative,  the  shift  operators  approach  can  be  used.  This  procedure  provides  a  lot  of 
different  master  formulas  depending  on  the  expression  taken  for  the  basic  integral.  We  will  just 
report  here  a  pair  of  ilustrative  examples. 

As  it  has  been  commented  previously  -see  sec  [3]-  shift  operators  allow  us  to  write  the  general 
integral  as: 


[  Xlm  I  Xvm'  ]  —  ^  ^ 


LML'M' 

k 


R)  J\ 


nLn'  L' 


(ill) 


k= 0 


where 


(R)  are  the  functions  defined  in  ec(p4h  and 


J\ 


nLn'  L' 


d_ 


ld_ 


d_ 

dC 


1  d  \L'  (1  d  \  L+L~k  .  i , 

'CWJ  \RdR  i 


(112) 


According  to  the  Fourier  transform  technique: 


“dt  cm) 


[  X  I  x'  ]  =  f(R ,  C,  0  =  [  Xoo  I  Too  ]  =  16  7T  VLJo  R  1/2  . 

Jo  (C  +  k2)  (C  +  k2) 

and  writing  the  derivative  operators  in  terms  of  Bessel  operators  the  application  of  the  operators 

is  trivial  and  gives: 


7 nLn'V  _  87T2  a/27T  (~1)L+L  k  ^  ^  „Ln(f\  JJri  ( t>\ 
Jk  ~  nL+L'-k+Ml  2^2^  *  ^1  H>  (S,  ) 
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J^L+L'-k+ 1/2 

l  l' 

kL+V-k- 1/2  jl+u 

-k+i,2(kR) 


dk 


(C2  +  k2)L+1+i  (C/2  +  k2)L'+l+i' 
where  E  <  i  <  n,  E  <i'<  n', and  the  functions  cfn( ()  are: 


^n(0  =  \lz 


[2  2l  n\  (. l  +  i)\  (2()2i~n  (-l)n+i 

7 r 


(114) 


(115) 


(n  —  i)\  (2 i  —  n)\ 

Notice  that  one-dimension  integrals  of  only  one  type  appear. 

The  semiinfinite  integrals  with  oscillatory  integrands  of  eq(  1 14  >  can  be  replaced  by  integrals 
with  finite  limits  and  nonoscillatory  integrands  by  means  of  Feynman  transform: 
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{p  +  p'  + 1)!  f1 


du 


up  { 1  -  u)p' 


(</2  +  k2)p+1  ((’2  +  k2Y+l  p\  pf\  J 0  (C2  u  +  C/2  (1  —  u)  +  k2)p+p'+2 

changing  the  order  of  the  integrals  and  taking  into  account  that: 
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(117) 


This  leads  to: 
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(20”  (20)”' 
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Jo  Jo 


where: 


im  = 


(_f)n+i  n|  (2C2  ^2)i 


(2i  —  n)!  ( n  —  i)\ 

The  analytical  solution  of  the  innermost  integrals  follows  immediately  from 


(119) 


dx  x2m+1  ku(ax )  = 
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2(a2/2)m+1  j- 


^(a2x2/2Yt 

/  J  .j  ^v+m+l— i \OLXj 


(120) 


and  the  second  integral  can  be  solved  with  numerical  methods. 

Eq(|l  18[)  allows  also  to  separate  the  long-  and  short-range  contributions.  The  former  is  attained 
by  bringing  the  upper  limit  of  the  innermost  integral  to  infinity,  and  the  latter  by  integrating  over 
1  <  t  <  oo.  Notice  that  the  long-range  contribution  must  coincide  with 


x 


(-1)''  a/vt  (l  +  l1  —  1/2)!  in  +  21  +  1)!  in'  +  2 V  +  1)! 
(Z  +  1/2)!  {V  +  1/2)!  (n+2l+ 2  Qn'+ 21' +2 

Im  I'm'  Zl+l'  (^) 
l+l' m"  pp2(l+l')+l 


(121) 


The  solution  in  case  of  equal  exponents  is  straightforward  by  taking  ( 
solving  the  pending  integral  with: 


('  in  eq(  114)  and 
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where  p  —  L  +  L' +  i  +  i'  + 1,1  —  L  +  L'  —  k  and  a  =  (. 

As  a  second  example,  we  will  consider  the  generalization  of 

}  (123) 

which  is  interesting  because,  in  this  case,  the  shift  operators  technique  allows  to  separate  the 
long-  and  short-range  contributions  as  well  as  to  obtain  simple  expressions  for  both. 

First,  we  write: 


[xW]  =  !6  vr2 
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Application  of  shift  operators  is  straightforward  since  the  variables  are  uncoupled,  and  gives: 
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Derivation  of  the  expression  for  the  short-range  term  is  a  bit  cumbersome.  This  term  is: 
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Bearing  in  mind  eq(|34]),  we  rewrite  this  term  as: 
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JJtLn  L  (short)  =  EEc"(o^'(0 
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(129) 


where  E  (^)  <  i  <  n,  E  <  i'  <  n'  and: 

=  (~i  r‘»!  (2Q2i- 

■ (n  -  >)!  (2 i  -  n)l  2‘  1  ' 

The  problem  is  thus  reduced  to  taking  the  three  pending  derivatives.  After  some  algebra,  it 
comes: 
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where  the  pending  integrals  are: 
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(132) 


7.1  Tests  on  the  accuracy  of  algorithms  for  Coulomb  integrals 

The  formal  developments  of  the  previous  sections  can  be  combined  in  many  ways,  leading 
to  a  huge  amount  of  possible  algorithms  for  the  calculation  integrals.  The  analysis  of  every 
alternative  implies  its  implementation  in  a  program  code  and  a  thorough  testing  of  the  results. 
As  a  consequence,  a  previous  selection  of  the  possible  candidates  is  mandatory. 

For  the  moment,  we  have  implemented  and  tested  only  three  algorithms  among  the  more  simple 
ones.  In  all  of  them,  the  final  Coulomb  integrals  are  obtained  as  linear  combinations  of  integrals 
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[  Xlm  I  Xl'm1  ]’  as  indicated  in  eq(  104).  The  algorithms  differ  in  the  way  the  [  Xlm  I  Xl'm >  ] 
integrals  are  computed. 

First,  a  code  in  double  precision  based  on  eq(]  1 1 2|)  combined  with  the  shift  operators  technique 
was  attempted  (file  coulomb _2  0 1 0_shiftop_D  .  f  90).  This  was  an  appealing  approach  be¬ 
cause  the  general  Coulomb  integrals  were  reduced  to  linear  combinations  of  overlap  integrals 
plus  one  Coulomb  integral  involving  functions  with  n  =  N  —  L  =  0.  Since  the  quantum 
numbers  appearing  in  the  overlap  integrals  are  twice  the  values  of  these  numbers  in  the  basis 
functions,  we  had  to  develop  a  new  algorithm  for  overlap  integrals  more  robust  than  any  other 
proposed  before.  Unfortunately,  after  the  new  algorithm  for  overlap  integrals  was  developed 
and  coded,  and  once  the  accuracy  in  the  overlap  was  guaranteed,  we  found  that,  for  high  quan¬ 
tum  numbers,  big  numerical  cancellations  occur  between  the  first  summand  in  eq(  1 12|)  and  the 
remaining  ones.  Because  of  it,  we  preferred  to  try  a  different  solution. 

In  the  second  attempt,  a  code  based  on  the  numerical  integration  of  eq(|110[)  was  prepared  (file 
coulomb  2010  intnum  D  .  f  90).  In  principle,  the  algorithm  seems  to  be  more  robust  than 
the  previous  one,  but  there  is  a  problem  with  the  dependence  of  the  integrand  with  the  values 
of  the  exponents  and  the  quantum  numbers.  This  dependence  is  rather  involved  and  makes  it 
difficult  to  design  a  good  choice  of  the  quadrature  points,  which  is  critical  in  the  final  result. 

In  view  of  this,  we  have  tried  a  third  solution,  which  is  closely  related  with  that  already  imple¬ 
mented  in  SMILES,  based  on  eq(|lQ8[).  The  algorithm  follows  closely  the  presciption  of  refll26l. 
but  the  calculation  of  the  auxiliary  functions  has  been  completely  redesigned.  In  this  way,  we 
have  been  able  to  extend  the  accuracy  of  the  calculations.  We  have  thus  prepared  three  codes 
corresponding  to  double  and  quadruple  precision  and  multiprecision [[271 .  The  latter  has  been 
taken  as  a  reliable  reference  for  all  the  remaining  codes. 

Figs|TT] to |~i~4| illustrate  the  quality  of  the  results  attained  with  the  double  and  quadruple  versions 
of  the  third  code  (with  the  multiprecision  values  taken  as  a  reference).  As  it  can  be  seen  in  figs 
[IT]  and  12,  for  high  quantum  numbers,  there  is  a  rather  broad  range  of  exponents,  defined  by 
the  value  of  the  lowest  ((',  Qr)  between  4  and  64,  in  which  there  is  a  significant  loss  of  accuracy. 
Using  quadruple  precision  in  this  range  is  sufficient  to  overcome  the  problem,  as  figs  [TT|  and  |T3| 
show. 

The  employment  of  quadruple  precision  implies  a  penalty  in  computational  time  by  a  factor 
between  10  and  100  with  respect  to  double  precision.  Since  the  algorithm  is  quite  fast,  the 
computational  cost  in  quadruple  still  remains  acceptable.  Nevertheless,  it  is  always  possible 
to  combine  both  double  and  quadruple  precision  in  the  same  code,  restricting  the  use  of  the 
quadruple  precision  to  those  cases  in  which  a  significant  loss  of  accuracy  has  been  detected. 
Finally,  it  should  be  noticed  that  the  formal  developments  above  reported  open  the  door  to 
many  other  possible  algorithms.  However,  at  this  moment,  we  consider  that  the  current  al¬ 
gorithm  based  on  translation,  with  quadruple  precision  in  the  prescripted  cases,  is  sufficiently 
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satisfactory  to  fulfill  the  current  requirements  in  molecular  calculations  with  STO. 
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Figure  11:  Accurate  decimal  figures  in  double  precision  for  Coulomb  integrals  with  N a  =  Nb  =  Nq  =  Njj  =  6, 
La  =  Lb  =  Lq  =  Ljj  =  5 
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Figure  12:  Accurate  significant  figures  in  double  precision  for  Coulomb  integrals  with  Na  =  Nb  =  Nc  =  No  = 

6,  La  =  Lb  =  Lq  =  Lb  =  5 
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Figure  13:  Accurate  decimal  figures  in  quadruple  precision  for  Coulomb  integrals  with  Na  =  Nb  =  Nc  = 
Nd  =6,  La  =  Lb  =  Lc  =  Lb  =  5 
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Figure  14:  Accurate  significant  figures  in  quadruple  precision  for  Coulomb  integrals  with  Na  =  Nb  =  Nc  = 
Nb  =6,  La  =  Lb  =  Lc  =  Lb  =  5 


8  Hybrid  integrals 


Hybrid  integrals  of  STOs  are  defined  as: 


'  nA  nB  |  nA  nA 

.  A-La  Ma  A- Lb  Mb  I  A.JJ  m'  X-L",  M", 


=  dr 


dr'  XZ  Ma  (Ca,  rA)  XnLBB  Mb  (Cb,  rB) 

1 


X  XLf  M'  (C A,  4)  XL»  m".  (C  4) 


r  —  r 


(133) 
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To  compute  the  hybrid  integrals,  the  one-center  distribution,  xjf  M,  {Q'a,  rA)  x'rf,  M„  (('" ,  rA)  is 
decomposed  as  a  linear  combination  of  one-center  distributions  -see  eq(|103 ).  In  this  way,  the 


general  hybrid  integral  of  eq(  133 )  can  be  expressed  as  a  linear  combination  of  integrals  like: 


[xlaaMa  xTbmb  |XA/J  =  dr  dr'xnLAAMA((A,r )  xTbMb^B^b) 


xL(C,r) 


(134) 

where  center  A  has  been  taken  as  the  coordinates  origin  and  Q  =  .  The  coefficients  of  the 

linear  combination  are  those  of  the  decomposition  of  products  of  regular  spherical  harmonics 
into  spherical  harmonics. 


The  problem  is  thus  reduced  to  the  efficient  computation  of  integrals  like  ( 134 ).  In  the  current 
approach,  these  integrals  are  computed  starting  from  some  basic  integrals  which  can  be  chosen 
as: 


ln. LM  =  J  dr  f  dr'xlAAMA((A,rA)  Xoo((B,rB)  XA/,(C,r')  (135) 

for  integrals  with  NB  even,  and 

JnLM  =  j  dr  f  dY'xnLM(U^A)  xIq{Cb^B)  ^  ^  Xx^C,*')  (136) 

for  Nb  odd,  and  translating  the  factor  |r  —  KAB\2E<'nB^2'>  (r -RB)  to  center  A.  In  a  lined-up 
axis  system  -with  the  Z  and  Z'  axes  coincident  and  the  (A",  Y)  axes  parallel  to  (A"',  Y')-  this 
translation  is  very  simple  and  all  the  integrals  with  M  ^  /i  vanish  after  integration  on  0. 

We  start  by  translating  the  factor  |r  —  Hab\2E^UbI2')  by: 

|r  -  RAB|2E(?lfl/2)  =  (r2  +  R\b  -  2rRAB  cos  9)E{nB/2]  (137) 

where  E(v)  stands  for  the  integer  part  of  v.  This  equation  yields  the  recurrence  relation: 

[n,  L ,  M;  vB  +  2,  0,  0  |  v.  A,  M ]  =  R2AB  [n,  L,  M;  vBl  0,  0  |  v,  A,  M } 

+  [n  +  2,  L,  M;  vB,  0,  0  |  u,  A,  M] 

~  Ul-M  +  1)  [n,L  +  l,M-,uB,0,0\u,\,M] 

+  (L  +  M)  [n  +  2,  L  —  1,  M;  uB,  0,  0  |  u,  A,  M]|-  (138) 

that  can  be  used  to  reach  the  integrals  for  uB  =  nB.  In  order  to  get  the  pertaining  values  of  LB 

and  Mb,  we  apply  the  translation  of  the  solid  regular  harmonic  (r  —  HB)  to  the  integrals 

[n,  L,  M;  nB,  0,  0  |  u,  A,  M],  by  means  of 
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(139) 


Mb 
ZT  B 

^ B 


(r- 


(-RAB)LB-k 


Mg 

Zk 


followed  by  the  decomposition  of  the  product  (r)  z^11  (r)  with 


4'( r)  zi? (r)  -  V  V  r21  ^+i,_2,(r)  (140) 

l 

where  |//|  <  l  <  E[(L  +  L')/ 2]  and,  in  the  lined-up  system,  the  sum  over  /i  contains  two  terms 
at  most. 

This  gives: 


[n,L,M-,nB,LB,MB  \u,X,fi]=  E  (  /f+ ijlM  ^  (~Rab)Lb  k 

k=\Mg\  ^  \  B\  / 

x  E  [rc  +  2/,  La  +  k  -  21,  ti-  nB,  0,0  I  u,X,n]  (141) 

i 

8.1  Calculation  of  the  basic  hybrid  integrals 

The  basic  integrals  can  be  written  as: 

Klm  =  j  drXnLAAMA^A,rA)  Xoo(C B,rB)  V^{(,r)  (142) 

and 

JnLM  =  I  drXnLAMA(CA,rA)  Xoo(C b,  Tfl)  V^(C,r)  (143) 

where  r)  is  the  potential  generated  by  the  one-electron  distribution  This 

potential  is  well  known  and  can  be  written  as: 


4  7 r 


f  (z/  +  2A  +  l)! 

1 

& 

1 

"i 

_ i 

)  £i/+ 2A+2  r2A+l 

7-| 

L  j=0  J'  J 

+ 


e"Cr 

C+1 


v-1 

£(c> 

3=0 


v\  {v  +  2A  +  1)! 
L7!  “  (2A+J  +  1)! 


Replacing 1 144| in p~42]  this  latter  can  be  written  as: 

ruXfl  _  \  ^  LM  Xfl 


tv^h  \  i^ivi  s\j n  a  nv 

1nLM  ~  2_^  aL+\-2l  0  ^ILX 


(144) 


(145) 


where: 
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Anv 

A ILX 


21  +  1  /("  +  2A  +  l)!ArI_x_1_1(aCB) 


X 


(2A  +  1)  C  l  C+2A+1 

2\  H 

J=0  J‘ 

zd  (1/  +  2A  +  1)! 


L.y!  U  •  2 A  •  1)! 

/if  (a,  3)  being  the  two-center  overlap  integrals: 


i/-i 


J2K+L+j+x~\(A  +  CXb) 

3=0 


(146) 


K(a,P)  =  /  dr  rn  zx  (r)  e 


Tb 


(147) 


For  the  J^lm  integrals,  an  equation  like  (]145[)  also  holds  replacing  hf  by  ///'  in  cq(]147[).  with: 


«r(a,«  =  y*^rf  (148) 

The  problem  is  thus  reduced  to  the  evaluation  of  the  overlap  integrals  hf(a,  3)  and  77f  (a,  3), 
which  can  be  accomphshed  with  any  of  the  procedures  previously  described  for  overlap  inte¬ 
grals.  In  particular,  the  current  codes  compute  these  integrals  by  means  of  the  STO  translation 
formulas: 


KM  3) 


4nR%l+1  e~(a+ ® 
21  +  1 


®i(P)  iin(oi,3)  + 


MP)  MM  3) 

(a  +  3)n+1 


(149) 


H?{a,p)  =  4t r  i^+i+2  e-(“+/3) 

1 


f  0 

'$l-l(P)  UnMP) 

M3)  il+l  nM  3Y 

21  +  1 

21-1 

21  +  3 

+ 


3  (a  +  3)n+1 


[MP)  h+ln  (a,  3)  -  <f>i-i(3)  hnMP)\ 


l  >  0  (150) 


MMP)  =  4n  Rn+B2  e-(a+V  U0n(a,  (3)-^  ilnMP) 


+ 


3  (a  +  3)n+1 


M3)  kln{a,3)  -  (  MP)  +  y  MP)  )  kon(a,3) 


MM 3)  =  yyy  [KMMP)  -  K+iMP)\ 


(151) 

(152) 


with 


<w)  =  ef  K,+1/M 

=  lF1(-l;-2l-,23) 


(153) 
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(154) 


fa{P)  =  (2/f3)l+1/2  (l  +  1/2)!  /;+1/2(/3) 

=  o-^i  ((  +  3/2;  f32 /A) 


iin(a,P)  =  (2/ p2)l+l/2  (. I  +  1/2)!  ea  [  dt  tn  e  at  (pt)l+1/2  Il+i/2(pt) 

Jo 

y-  (/32/ 4)J  n  +  2/  +  2j  +  3;  a) 
j!  (/  +  3/2) j  (n  +  21  +  2j  +  2) 


hn(ot,0) 


,a+0  {a  +  /9)n+1  f!  2;+1/2 

(2/)!  7tV2 


Cn+j  (O!  +  /5) 


dt  tn  e~at  (/3t)/+1/2  A)+1/2(/3t) 

H)j  (w  + j)! 

j!  (-20; 


(156) 


In  these  equations,  (a)j  denotes  the  corresponding  Pochhammer  symbol,  Kv(z)  are  the  Mac¬ 
donald  functions,  Iv(z)  the  corresponding  Bessel  functions,  and  en(z)  stands  for  the  truncated 
exponential: 


71  3 

Zn(z)  =  Y.~\  (157) 

1=0  J ' 

The  set  of  functions  is  computed  by  recursion  in  a  fully  stable  way: 

»WW=«,W+(2|-.1f(2,  +  1)».-.(W  (158) 

starting  from  <T>0(/?)  =  1  and  $i(/?)  =  1  +  p. 

The  set  of  fa  functions  can  be  computed  also  by  recursion  by  means  of: 


fa-M  =  MP)  + 


(J2 


fa+l(p) 


(159) 


(2 1  +  1)  (2 1  +  3) 

since  the  relation  is  stable  for  backwards  recursion,  it  must  be  started  from  famax(P)  and  famax-i(P), 
which  are  computed  by  means  of  eq(  154).  To  improve  performance,  Miller’s  algorithm[|28l  is 
used. 

For  the  calculation  of  the  kin(a,  P),  the  first  row  is  computed  by: 


the  second  row,  by: 


k0n(a,P)  =  n\  en(a  +  P) 


(160) 
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(161) 


kln(a,/3)  =  k0n(a,p)  H - k0n+1(a,P) 

a  +  p 


and  the  remaining  ones,  by: 

ki+ln(a,p)  =  kin(a,/3 )  + 


P2 


k, 


l — Iti — | — 2 


[a,p) 


(162) 


(a  +  /?)2  (2Z  +  1)  (21  —  1) 

Finally,  in  the  calculation  of  the  iin(a,  P),  three  cases  are  distinguished: 

Case  1:  P  <  1 

The  elements  iimaxn(ot,  P)  with  n  <  —  1  are  computed  by  means  of  the  second  expression  in 
eq(  155|).  Although  the  sum  is  infinite,  it  converges  very  quickly  for  these  values  of  6.  The 
remaining  elements  are  computed  by  means  of: 


iin(a,P)  = 


21  +  3 


i+i(P )  +  a  P+ in-i(a,  p)-(n-  1)  ii+ in_2(a,  p)\ 


(163) 


P  n+ 1  j  /  -l) 


(21  +  3)  (n  +  21  +  3) 


(a2  -  P2)  ii+ln+1(a,p)  -  (n+  1)  a  ii+ln(a,P) 


+  (a  +  21  +  3)  pi+i(P)  + 


P2 


21  +  5 


Pl+2(P) 


(164) 


iin(a,P)  = 


{(a2  -  P2)  a  ii+ln+i(a,P ) 


(21  +  3)  (n  +  21  +  3)  (n  +  21  +  2) 

—  [(n  +  1)  a2  +  (n  +  2Z  +  3)  /32]  ii+in(oc,  P)  +  cx.~  pi+i(P) 

+  (21  +  3)  (a  +  n  +  21  +  3)  4>i(P) } 


(165) 


Case  2:  (3  >  a  +  8 

The  elements  ii-2imax(c<,  P)  and  11-21^+1(01,  P)  are  computed  by  eq(|155|).  The  remaining 
elements  with  l  =  lmax  are  computed  by: 


P  n+2  (ot,  P) 


a 2  —  P2 


~(n  +  1)  (n  +  21  +  2)  iin(a,  P) 


+  2  a  (n  +  21  +  2)  iin+\(ot ,  P)  +  (n  +  1  —  a)  4>i(P) 


P2 


<j>i+1(P)(l66) 


21  +  3 

Case  3:  1  <  (3  <  a  +  8 

The  elements  %in((y,1  P)  with  ti  arnjn  and  1  a.m.ax  ^ mint  <ffld  /  'arn(,x  1 


are  computed  by  eq(  155 ).  The  elements  with  n  =  nmin  +  1  and  l  =  nmax  —  nmin  —  1,  and 
l  =  nmax  -  nmin  -  2  are  computed  by: 


iin+i(a,P)  =  - 
a 


(n  +  21  +  2)  iin(&,  P)  + 


P2 


21  +  3 


ii+ln(a,P)  -  4>i(P) 


(167) 
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and  ([T64])  respectively. 

The  remaining  elements  with  n  =  nmm  and  n  =  nmm  +  1  are  computed  by: 


iin(a,(3)  = 


(n  +  21  +  2)  ( n  +  21  +  3) 

(«2  -  / 3 2)  P2 


(■ n  +  21  +  7/2)  p2 
1  +  3/2 


a 


ii+in(a,P) 


(21  +  3)  (21  +  5) 

and  the  remaining  elements,  by: 


h+2 n(os  P)  +  (a  +  n  +  21  +  3)  4>i(P)  + 


P2 


21  +  3 


<f>i+i(P)  (168) 


il-ln+2(ot,P)  =  iin(a,P)  + 


p2 


(21  +  3)  (21  +  1) 


ii+in(a,P) 


(169) 


Finally,  the  hypergeometrics  iFi(l;  n;  a)  =  Fn(a )  appearing  in  eq(  155 )  are  computed  by  back¬ 
wards  recursion: 


Fn(oc )  -14 — -  Fn+i(a) 
n 

starting  with  the  explicit  definition 


(170) 


OO  A 

_  /  x  v — ^  Or 

Fn(<y)  =  Y,~  (171^ 

J=0  J ' 

for  an  index  n  such  that  the  series  converges  quickly. 

8.2  Tests  on  the  accuracy  of  algorithms  for  hybrid  integrals 

The  above  exposed  algorithms  for  hybrid  integrals  have  been  implemented  in  FORTRAN  at 
three  different  levels  of  accuracy:  double  precision,  quadruple  precision  and  multiprecision Ii27l . 
The  multiprecision  version  (with  a  working  precision  of  65  decimal  digits)  has  been  used  as  a 
reference  for  testing  the  accuracy  of  the  other  two. 

Figs  [13]  to  [20]  illustrate  the  accuracy  attained  in  the  hybrid  integrals  with  the  current  algorithm 
for  N  =  6  and  L  =  5  for  the  four  functions.  As  it  can  be  seen,  there  is  a  wide  range  of  values 
of  the  exponents  for  which  the  accuracy  loss  in  double  precision  is  dramatic,  and  the  algorithm 
seems  not  no  be  suitable  in  DR  When  quadruple  precision  is  used,  the  range  becomes  narrower 
but  it  still  remains  a  range  -corresponding  to  values  of  the  exponents  that  do  not  appear  in  usual 
basis  sets-  for  which  all  the  figures  are  lost  even  in  quadruple  precision.  The  consequence 
is  rather  obvious,  if  exponents  within  this  range  will  be  used,  algorithms  based  on  a  different 
formalism  must  be  investigated.  For  the  moment,  the  only  available  solution  is  based  in  the 
Gaussian  expansions  of  the  STO,  which  can  be  used  to  cover  these  pathological  cases. 
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Hybrid  integrals:  exA=exB;  exA'=exB' 

Accurate  decimal  places  NA=NA'=6  LA=LA'=5  NB=NB'=6  LB  =  LB'  =  5 
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Figure  15:  Hybrid  integrals  with  C^a  =  C b  and  Q'A  —  C  B  in  double  precision 


Hybrid  integrals:  exA'=exA";  exA'=exB+5 
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Figure  16:  Hybrid  integrals  with  C,a  =  (b  and  C'b  =  (a  +  5  in  double  precision 
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Hybrid  integrals:  exA'=exA";  exB=exA+5 
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Figure  17:  Hybrid  integrals  with  Cb  =  Ca  +  5  and  C,'A  —  Cb  +  5  in  double  precision 


Hybrid  integrals:  exA=exB;  exA'=exB' 

Accurate  decimal  places  NA=NA'=6  LA=LA'=5  NB=NB'=6  LB  =  LB'  =  5 
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Figure  18:  Hybrid  integrals  with  Qa  =  and  Q'A  =  ('B  in  quadruple  precision 
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Hybrid  integrals:  exA'=exA";  exA'=exB+5 
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Figure  19:  Hybrid  integrals  with  Qa  =  C b  and  =  Q'A  +  5  in  quadruple  precision 


Hybrid  integrals:  exA'=exA";  exB=exA+5 
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Figure  20:  Hybrid  integrals  with  Qj  =  <^4  +  5  and  ('A  =  ('B  +  5  in  quadruple  precision 
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9  Exchange  integrals 

Exchange  integrals  of  STOs  are  defined  as: 


J  dr  J  dr' xIaaMa(Ca,  ta)  xTbmb  (Cb^b) 


(172) 


Although  the  translation  method  previously  exposed  can  be  also  used  for  these  integrals,  carry¬ 
ing  out  a  double  translation,  in  this  case  the  treatment  in  terms  of  ellipsoidal  coordinates  seems 
to  be  better  suited,  and  the  algorithm  currently  implemented  is  based  on  this  approach.  To  our 
knowledge,  the  first  satisfactory  solution  to  these  integrals  using  ellisoidal  coordinates  was  re¬ 
ported  by  Ruedenberg  ll29ll.  who  proposed  a  factorization  of  the  integral  as  a  product  of  two 
functions,  each  associated  to  a  charge  distribution.  The  procedure  involves  a  single  numerical 
integration  and,  unfortunately,  fails  to  give  even  middle  accuracy  for  high  quantum  numbers  and 
small  exponents.  The  presence  of  the  numerical  integral,  though  simple,  makes  it  difficult  to 
cover  all  the  combinations  of  exponents  and  quantum  numbers  in  a  simple  way.  An  alternative 
was  proposed  by  Maslen  and  Trefry  lf30ll  which  leads  to  some  improvement  in  the  performance, 
but  with  quite  a  significant  downgrading  in  the  accuracy  when  increasing  quantum  numbers. 
We  propose  a  better  procedure  based  on  ellipsoidal  coordinates  and  the  extensive  usage  of 
recurrence  relations  [t3TI  is  summarized  below. 

In  this  algorithm,  the  exchange  integral  are  written  as: 


n-A+LA+nB  nA+^A+nS 
+LB-\m\  JrL'B  —  \m\ 


X 


k= 0 


^  jUaLaIMaI^bLeIMbI  J^A^aI^aI^'b^bI-^bI  (ryf\ 



(173) 


where 


M+  =  sgn (MaMb)  x  (\Ma\  +  \Mb\) 
M_  =  sgn (MaMb)  x  |  \Ma\  -  | Mg | 


(174) 
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w'r\k,k'AA')  = 

/»oo 

J  it;  (f 2  -  1)”/2  Q?(0  [£*  e-«  K£ (P',0 

+ 

f  K?k(P,0] 

(175) 

= 

J  d£  (£/2  -  l)m/2  P™{xi')  e~x(: 

(176) 

jnLM,n'L'M' 
Jlkm  1 

n+n'+L 
+L'— \m\ 

:*)  =  «-  y.  w 

(177) 

3=  0 


where  sm  and  L  A/ '  are  numerical  coefficients  whose  calculation  is  summarized  in  the 

appendix  B,  and 

irM  =  j  <*1  (i  -  ’?2r/2  prw  <f  wo 

Finally: 


7  —  (C a  ~  (b)  Rab/2  (3  —  ((a  +  C b)  Rab/2 
T,  =  (Ca-S'b)Rab/ 2  F  =  (C'a  +  C'b)Rab/ 2 


It  can  be  seen  that  both  the  functions  J, 


ikm'n  L  M  (z)  and  W\m\ki  k!\ ( 3 ,  (3')  play  an  essential 
role  and  the  efficiency  of  the  the  algorithm  is  determined  by  their  calculation. 

According  to  eq(  177 ),  the  calculation  of  the  L  M  (z)  requires  the  L  A/ '  coef- 


l  km 

ficients  and  the  i)™(z).  Since  the  L  A/ '  do  not  depend  on  the  exponents  of  the  STO, 

they  can  be  computed  once,  following  the  prescription  of  appendix,  and  stored.  The  0*1  are 


Lip ' 


closely  related  to  the  4>i(z)  functions  defined  in  eq(  154),  and  their  efficient  computation  follows 
the  same  scheme  as  described  for  those  functions. 

Nevertheless,  it  should  be  noticed  that  the  L  M  (z)  depend  on  the  parameters  of  the  distri¬ 
butions  (pairs  of  STO)  whereas  the  L  A/ '  depend  on  pairs  of  distributions  (four  STO). 

As  a  consequence,  the  efficiency  of  the  overall  algorithm  is  determined  by  the  computation  of 
these  functions.  The  algorithm  for  their  efficient  calculation  is  described  below. 


9.1  Calculation  of  the  functions  wjm^  (. k ,  k'\  / 3 ,  f3') 

The  recurrence  relation  to  increase  m: 


54 


wr+i(k,k'-,p,p') 


(/  —  m) 


{l  +  1  —  771 )2 


wl"l  (fc,  fc';  ft/?) 


21  +  1 

(l  +  m  +  1)  W\m\k  +  l,k'  +  l-f3,!3') 


(179) 


reduces  the  problem  to  evaluating  the  set  of  functions  with  m  =  0.  To  compute  them,  it  is  better 
to  work  with  the  auxiliary  functions: 


Obviously: 


mkk'  = 
WIV  — 


+  /  <RQi'RK‘'e-«  /  <%' fttf)  e*  e-K 


(180) 


W,W  (fe.fc';  ft  /?')=<’  (181) 

From  the  recurrence  relation  of  the  Legendre  polynomials  -see  0  eq  8.733.2-  the  following 
recurrence  relations  for  the  can  be  derived: 


with 


(2f  + 1)  w++i  =  (r  + 1)  t«j?r+1 

(2i+i)mf+i‘'  =  (i+ix;1, 
=  K.,?'-1  +  7  [iW-iW  +  /?) 
w+  =  w^lk'  +  1  [KwiiP  +  P) 


l'  u’ifi-1  ('>0  (182) 

l  wkk'u,  l  >  0  (183) 

(A/  -  1)  wf/”1  -  w**'"2]  (184) 

(fc-l)Wf0-lfc'-<-2fc']  (185) 


Kln(s)=Kl(s,  oo)  (186) 

From  these  equations,  it  follows  that  the  set  of  functions  u^fc'  can  be  computed  in  a  numerically 
stable  way  from  the  functions  ww  =  u;))0 .  For  these  functions,  the  following  relations  hold: 


ww 


P 

2Z'  +  1 


[ww+i  ~  U’ll'-l 


+  Uu> 


(187) 
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where 


ww 


(3 

21  +  1 


\wi+H>  ~  Wi-n' 


+  Ui'i 


(188) 


/•LXJ 

uw  =  uu,(s)  =  /  e-*  [Q,(£)  Pt1^)  -  m)  Q{i'\xi)} 


(189) 


with 


p/-1’®  =  (V  p,(  a  =  ^  [p,+i(?)  -  p, -,«)] 

Ql1)(0  =  J‘ V Q|(0  =  [<2‘+1K) -*-!«)]  ,>0 


(190) 


(191) 


The  recursions  ( 187 )  and  ( 188 )  are  stable  neither  for  ascending  nor  for  descending.  Because  of 
it,  the  bisection  algorithm  reported  in  Il32ll  has  been  used,  starting  from  w 00,  Woimax,  wimax o  and 
wimaximax.  The  first  three  are  computed  by: 


moot, 3,0')  =  ^{ed-d' T(O,20) +ed,-dT(0, 20)  ^ed+d'T[0, 2(0  +  0')] 


+  e 


-U3+0') 


C  +  ln 


W 
(3  +  P' 


(192) 


C  being  Euler  constant  (0.577  215  664  ...); 


wio(f3,f3')  = 


with  l  =  lmax,  and 


1 

P 


e  Li((3)  +  vi- i  - 


eW 


i- i 


£ 


(200 

k  +  1 


6?(f3)T{-k,2(P  +  p')\ 


0i(/3)  M/3',/3  +  f3') 


wqi({3,  (3')  =  wio(f3',  (3) 


(193) 

(194) 


T lmax  (/^) 


rv 


di  e“* 


lmax  ( lmax  + 1)  L  ((  max 


l  - 


=  E 


((+7)! 


E  (7  —  j)\  j\  (2/3)7 


_ 

1)  (lmax  +  1) 


(195) 

(196) 


56 


(197) 


OO  /  (~y  \k 

vi{x,  y)  =  ey  T(y~kl  2y } 


'OO 


dt  e~l  t-71-1 


(198) 


For  f3/j3'  <  5,  the  following  equation  is  used  instead  of  ( 197 ): 


-  e^^fr{-k,2y) 

k= 0 


(199) 


The  value  of  lmax  is  fixed  as  a  power  of  2  large  enough  to  allow  to  compute  the  wimaximax  by  the 
first  term  of  its  asymptotic  series: 


e~<j3+P’) 


(200) 


wu~  (p+wni+i)* 


9.2  Tests  on  the  accuracy  of  algorithms  for  exchange  integrals 

The  above  exposed  algorithm  for  exchange  integrals  have  been  implemented  in  FORTRAN  at 
three  different  levels  of  accuracy:  double  precision,  quadruple  precision  and  multiprecision IG71 . 
The  multiprecision  version  (with  a  working  precision  of  65  decimal  digits)  has  been  used  as  a 
reference  for  testing  the  accuracy  of  the  other  two. 

Figs  [2T| to  [24] illustrate  the  accuracy  attained  in  the  exchange  integrals.  As  it  can  be  seen,  there 
is  a  region  of  intermediate  values  of  the  exponents  in  which  some  loss  of  accuracy  is  found, 
mainly  when  the  exponents  of  the  functions  of  a  given  distribution  are  different.  In  fact,  it  is  to 
be  expected  that  higher  the  differences  between  the  exponents  imply  a  slower  convergence  on 
the  series  on  l  of  eq(|173|).  Nevertheless,  the  accuracy  loss  in  the  cases  tested  seems  to  be  not 
too  dramatic  and,  if  higher  accuracy  is  required,  the  quadruple  precision  does  the  job  in  a  very 
satisfactory  way,  as  it  can  be  appreciated  in  figs  10  to  12.  As  a  consequence,  most  integrals  can 
be  computed  in  double  precision  and  the  use  of  quadruple  precision  can  be  restricted  to  those 
cases  in  which  a  significant  loss  of  accuracy  has  been  detected. 

Finally,  it  should  be  mentioned  that  algorithms  for  exchange  integrals  based  on  the  translation 
of  STO  may  be  an  interesting  alternative  in  cases  involving  charge  distributions  with  a  high 
asymmetry,  i.e  with  one  exponent  much  larger  than  the  other. 
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Exchange  integrals:  exA=exB;  exA'=exB' 

Accurate  decimal  places  NA=NA'=6  LA=LA'=5  NB=NB'=6  LB  =  LB'  =  5 


Double  precision 
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Figure  21:  Exchange  integrals  with  (a  =  C b  and  C,'A  =  CB  in  double  precision 


Exchange  integrals:  exB=exA;  exB'=exA'+5 

Accurate  decimal  places  NA=NA'=6  LA=LA'=5  NB=NB'=6  LB  =  LB'  =  5 

Double  precision 

exA'=exB1 0.0025 
exA 

0.005 

0.01 

0.03 

0.06 

0.12 

0.25 

0.5 

1 

2 

4 

8 

16 

32 

64 

0.005 

32 

31 

31 

31 

30 

30 

30 

29 

29 

29 

29 

29 

29 

33 

42 

0.01 

30 

29 

29 

29 

28 

28 

28 

27 

27 

27 

27 

27 

27 

30 

41 

0.03 

27 

26 

26 

26 

25 

25 

25 

24 

24 

24 

24 

24 

24 

27 

38 

0.06 

25 

24 

24 

24 

23 

23 

23 

23 

22 

22 

22 

22 

22 

26 

36 

0.12 

23 

23 

22 

22 

22 

21 

21 

21 

20 

20 

20 

20 

20 

24 

34 

0.25 

21 

21 

20 

20 

20 

19 

19 

19 

19 

18 

18 

19 

19 

22 

32 

0.5 

19 

19 

19 

18 

18 

18 

17 

17 

17 

17 

17 

17 

16 

20 

30 

1 

18 

18 

17 

17 

17 

16 

16 

16 

16 

15 

15 

15 

15 

19 

28 

2 

17 

17 

16 

16 

16 

16 

15 

15 

15 

14 

14 

13 

13 

17 

27 

4 

17 

16 

16 

16 

15 

15 

15 

14 

14 

14 

14 

12 

12 

15 

25 

8 

17 

16 

16 

16 

15 

15 

14 

14 

14 

14 

12 

11 

11 

14 

25 

16 

17 

16 

16 

16 

15 

15 

15 

14 

14 

14 

13 

12 

12 

15 

25 

32 

20 

20 

20 

19 

19 

19 

19 

18 

18 

18 

17 

15 

15 

19 

29 

64 

31 

31 

30 

30 

30 

29 

29 

29 

29 

29 

27 

26 

26 

29 

39 

Figure  22:  Exchange  integrals  with  (A  =  C b  and  Q'B  =  C,'A  +  5  in  double  precision 
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Exchange  integrals:  exB=exA+5;  exA'=exB'+5 
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Figure  23:  Exchange  integrals  with  C-B  =  Ca  +  5  and  Ca  =  Cb  +  5  in  double  precision 


Exchange  integrals:  exA=exB;  exA'=exB' 
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Figure  24:  Exchange  integrals  with  Qa  =  (b  and  C 'a  =  Cb  m  quadruple  precision 
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Exchange  integrals:  exB=exA+5;  exA'=exB' 
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Figure  25:  Exchange  integrals  with  £4  =  Qj  and  C,'B  =  (,'a  +  5  in  quadruple  precision 


Exchange  integrals:  exB=exA+5;  exA'=exB'+5 
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Figure  26:  Exchange  integrals  with  Cb  =  Ca  +  5  and  (,'A  =  C,'B  +  5  in  quadruple  precision 
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10  Large  STO-nG  expansions  for  three-  and  four-center  two-electron  in¬ 
tegrals 

Up  to  date,  the  efforts  carried  out  for  the  direct  calculation  of  three-  and  four-center  two-electron 
integrals  involving  STO  have  not  yield  sufficiently  efficient  algorithms  yet.  Although  it  is  to  be 
expected  that  these  efforts  will  prove  fruitful  in  the  near  future  and  that  such  efficient  algorithms 
will  be  designed,  the  only  currently  available  algorithms  for  computing  these  two-electron  in¬ 
tegrals  in  the  general  case  are  based  on  the  use  of  large  Gaussian  expansions  of  the  STO  -the 
so-called  STO-nG  expansions. 

Using  these  expansions,  each  integral  with  STO  becomes  a  linear  combination  of  integrals 
involving  Gaussian  type  orbitals  (GTO),  which  can  be  computed  by  very  efficient  algorithms. 
In  order  to  attain  an  accuracy  in  the  integrals  suitable  for  molecular  calculations,  the  number 
of  Gaussians  per  STO  can  be  moderately  large  (between  10  to  15  in  HF  calculations,  larger  for 
Cl  calculations)  and  usual  algorithms  implemented  in  standard  computational  packages  cannot 
be  very  efficient  in  this  case.  In  1996  Prof.  K.  Ishida  developed  a  new  algorithm[f33l.  and  has 
improved  it  later  on.  This  algorithm  has  proved  to  be  very  efficient  when  using  large  STO-nG 
expansions,  and  it  is  currently  implemented  in  SMILESII21. 

In  the  final  phase  of  the  current  project,  new  STO-nG  expansions  have  been  attained  to  cover  a 
range  of  quantum  numbers:  0  <  L  <  6  and  L  <  N  <9  except  the  Os)  that  is  sufficient  for  most 
cases.  Attaining  long  expansions  is  not  an  easy  task,  mainly  for  high  quantum  numbers,  and  it 
requires  special  techniques  for  minimization  of  residual  funtions  depending  on  a  large  number 
of  nonlinear  (exponential)  parameters.  In  our  group  we  have  developed  one  of  such  techniques, 
which  is  summarized  in  appendix  C. 

The  expansions  obtained  in  this  project  have  been  added  to  the  file  with  the  previously  existing 
ones,  and  the  new  file  is  included  in  the  supplementary  material  of  this  report,  together  with  a 
standalone  program  for  computing  two-electron  integrals  with  these  expansions. 

11  Conclusions 

Different  algorithms  for  integrals  appearing  in  molecular  integrals  with  Slater  type  orbitals 
(STO)  have  been  coded  and  their  performance  has  been  thoroughly  tested  regarding  their  accu¬ 
racy  and  computational  cost. 

In  a  first  step,  algorithms  for  the  calculation  of  two-center  one-electron  integrals  have  been 
formulated,  coded  and  analyzed.  The  main  conclusions  of  this  part  were  that: 

•  The  algorithms  already  currently  available  in  SMILES  for  these  integrals,  based  on  recur¬ 
rence  relations,  cannot  provide  sufficient  accuracy  for  large  N  and  L  quantum  numbers. 
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•  The  new  algorithm  using  ellipsoidal  coordinates  enables  to  compute  these  integrals  with 
sufficient  accuracy  for  molecular  calculations  in  a  range  of  values  of  N  and  L  considerably 
larger  than  the  previous  one. 

•  An  alternative  algorithm  based  on  the  shift  operators  technique  has  been  also  developed 
and  tested  to  be  used  in  some  algorithms  for  the  computation  of  Coulomb  integrals. 

In  a  second  step,  the  calculation  of  the  two-center  coulomb  integrals  has  been  dealt  with.  Since 
these  integrals  play  a  very  important  role  in  some  formulations  of  molecular  structure  calcu¬ 
lation,  a  thorough  analysis  of  the  mathematical  aspects  of  the  problem  was  carried  out.  The 
formulation  in  terms  of  the  shift-operators  technique,  Fourier  transform  and  translation  method 
was  carefully  revisited  and  different  combinations  of  these  techniques  were  explored.  Finally, 
three  different  algorithms  were  selected  for  implementation  and  testing:  one  based  on  shift  op¬ 
erators,  and  two  different  algorithms  using  the  translation  of  the  STO.  The  main  conclusions  of 
this  part  were: 

•  The  extension  to  higher  quantum  numbers  of  the  algorithms  already  available  in  SMILES 
for  theseintegrals,  based  on  the  translation  of  STO,  implies  further  parametrization  of  aux¬ 
iliary  functions,  which  is  a  cumbersome  process. 

•  The  new  algorithm  using  shift  operators  leads  to  serious  numerical  cancellations  for  high 
quantum  numbers  which  cause  an  important  loss  of  accuracy. 

•  The  new  algorithm  using  translation  methods  accompanied  by  a  one-dimension  numerical 
integration  seems  to  be  more  robust  and  is  quite  fast,  but  the  dependence  of  the  integrand 
with  both  the  exponents  and  quantum  numbers  of  the  STO  makes  it  difficult  to  implement 
an  efficient  quadrature  scheme  for  the  general  case. 

•  The  algorithm  using  translation  methods  with  analytical  integration  seems  to  be  the  best 
option  in  this  moment. 

In  the  third  step,  the  hybrid  integrals  were  examined  and  algorithms  based  on  the  translation 
method  were  implemented,  coded  and  tested.  The  main  conclusions  of  this  part  were: 

•  The  algorithm  used  in  SMILES  for  hybrid  integrals,  based  on  the  translation  of  STO,  lead 
to  a  serious  loss  of  accuracy  for  high  quantum  numbers  in  the  current  implementation. 

•  The  new  algorithm  for  hybrid  integrals,  based  also  in  the  translation  method  but  with 
different  procedures  for  the  computation  of  auxiliary  functions,  allows  to  extend  the  range 
of  applicability  to  higher  quantum  numbers,  but  quadruple  precision  is  mandatory  in  a 
broad  range  of  values  of  the  exponents,  this  being  insufficient  in  some  combinations  of 
large  and  small  values  of  the  exponents. 
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•  Algorithms  based  on  the  usage  of  large  STO-nG  expansions  of  the  STO  (see  below  in  this 
report)  can  be  a  provisional  alternative  for  these  pathological  cases,  but  a  further  study  of 
alternative  techniques  such  as  those  based  on  ellipsoidal  coordinates,  seems  to  be  desir¬ 
able. 

The  third  step  was  devoted  to  the  computation  of  exchange  integrals  with  algorithms  based  on 
ellipsoidal  coordinates.  The  main  conclusion  of  the  analysis  were: 

•  The  algorithm  used  in  SMILES  for  exchange  integrals,  based  on  the  translation  of  STO, 
can  be  confidently  used  for  moderately  high  quantum  numbers,  but  changes  must  be  done 
in  the  codes. 

•  The  new  codes  for  these  integrals  allow  to  extend  the  range  of  applicability  to  higher  quan¬ 
tum  numbers,  the  bottleneck  being  now  the  size  of  some  intermediate  auxiliary  matrices, 
which  hinders  the  usage  of  multiprecision.  Nevertheless,  this  is  not  a  serious  drawback 
since  multiprecision  is  only  used  as  a  testing  reference. 

•  A  further  effort  to  redesign  the  implementation  avoiding  the  use  of  such  large  intermediate 
matrices  may  be  convenient. 

Finally,  as  scheduled  in  the  proposal,  new  STO-nG  expansions  were  developed  to  fully  cover 
the  range  of  quantum  numbers:  0  <  L  <  6  and  L  <  N  <  9  (except  the  Os).  These  expansions 
are  used  for  computing  three-  and  four-center  molecular  integrals  with  STO  and,  as  commented 
above,  could  be  useful  to  cover  the  pathological  cases  of  hybrid  integrals  until  a  better  alternative 
is  available. 

Instructions  for  installing  and  running  the  programs  corresponding  to  the  different  algorithms 
accompanying  this  report  are  given  in  appendix  D. 


Appendix  A.  The  functions  Aj(j3)  and  Bj{y ) 


The  functions  Aj(/3)  are  closely  related  to  the  Euler  Incomplete  Gamma  functions  of  second 
type: 


r°°  p 

AM  =  J1  ?  e~“  =  ^  etf) 

where  ej(/3)  is  the  expansion  of  the  exponential  function  truncated  at  jth  order. 
They  can  be  computed  in  a  fully  stable  way  by  ascending  recursion: 
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The  calculation  of  the  Bj(u)  functions  is  a  bit  more  difficult.  Their  main  properties  are  summa¬ 
rized  below. 

Definition: 


Bj(u)  =  j  drj  rf  e  vr} 

Power  expansion , 
for  j  even  (j  =  2 p): 


(203) 


B2p{y)  = 


,2  i 


^(p+l/2  +  i)(2i)!  P+1/2 


i F2(p  +  1/2;  1/2,  p  +  3/2;  i/2/4)  (204) 


for  j  odd  (j  =  2p  +  1): 
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Expansion  in  Bessel  functions: 

for  j  =  2 p  the  following  identity  can  be  used: 
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+  ’  2s  ip  —  n)\  n\  1  j 
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and  therefore  (see  0  eq.  8.431): 
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For  j  =  2p  +  1,  combining  the  relation: 
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B2p+ \{v)  =  -^B2p(v) 


with  0  eq.  8.486.5,  it  follows: 
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Expansion  in  Gamma  functions. 


Direct  integration  of  eq(  203 )  yields: 
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(210) 


Bj{u)  =  2ev  /  dt  (2t  -  l)j  e~2vt  = 
Jo 

Recurrence  relations. 

Integration  by  parts  gives: 
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This  relation,  applied  twice,  gives: 
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Iterated  application  of  this  process  yields: 
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'I'he  algoritm  for  Bj{v)  functions. 

The  algorithm  followed  to  compute  the  Bj(y)  functions  consists  of  the  use  of  upwards  recursion 
from  j  —  0  to  j  —  E[y]  and  downwards  recursion  from  j  =  jmax  to  j  =  E[v].  Upwards 
recursion  starts  with: 


Bn(u)  = 


2  sinh(y) 


v 


B0{Q)  =  2 


Downwards  recursion  starts  with  BJrnax  (y)  computed  by  means  of  eq(  213)  with  jr 


(216) 
taken  as 


the  lowest  suitable  even  number  to  prevent  cancellation  errors.  Convergence  of  the  series  is 
checked  and  a  warning  message  will  be  issued  if  the  prescribed  convergence  (l.d-35  in  quadru¬ 
ple  precision)  is  not  achieved  with  the  highest  allowed  number  of  terms  (100).  We  have  never 
observed  that  message  throughout  our  tests,  i.e.  the  maximum  available  number  of  terms  seems 
to  be  sufficient  for  reaching  convergence  in  all  cases. 
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The  stability  of  the  recursion  scheme  was  thoroughly  tested  with  Mathematica,  working  in 
extended  precision,  and  no  loss  of  accuracy  was  observed.  The  algorithm  seems  to  be  fully 
stable. 


Appendix  B.  Calculation  of  the  coefficients  trn  -  ’ 

The  decomposition  of  a  two-center  charge  distribution  in  ellipsoidal  coordinates  gives: 


where 


with 


n+n'+L+L' 


’•a  rB  xIm^a,  0  xi'M’iyB, o  =  (f'f'j 

x  e-v  R"/2  Yi  sAW  [<e2  -  1)  (1  -  i?2)] 


e 

2\1  M/2 


-£  (C+CO  Rab/2 


X 


n+n'+L  n+n'+L 
—  \fi\  +L1  —  |/z| 

E  E  « 

r=0  s=0 


nLIMI.n'L'IM'l  s 


(217) 


cos  [if  fi  >  0 
sin  |/r|0  [i  <  0 


(218) 


-V  3^(0)  —  ^Ar(0)  $M'(0)  —  %+  ^M+(0)  +  Sm_  (219) 


{1  if  M  =  0  or  M'  =  0 

—1/2  if  sgn(M)  <  0  and  sgn(M')  <  0  (220) 

1/2  otherwise 

{0  if  M  =  0  or  M'  =  0 

1/2  if  sgn(M)  =  sgn(M')  (221) 

sgn(M)  sgn(|M'|  -  \M\)  (1  -  5M-m>, o)/2  otherwise 

M+  =  (|M|  +  \M'\)  sgn (M)  sgn(M')  (222) 

M_  =  ||M|  -  \M'\  \  sgn (M)  sgn(M')  (223) 

where  sgn (m)  stands  for  the  algebraic  sign  of  rn  times  1. 

The  L  M  coefficients  can  be  expressed  in  terms  of  and  /i"  A!'U  defined  by: 
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which  can  be  computed  by  recursion: 
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The  final  expression  of  the  q^lm A V M'  being: 
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Appendix  C.  Approximation  of  functions  with  basis  sets  containing  a  big 
number  of  nonlinear  parameters 

Let  us  consider  the  approximation  of  a  function  /( r)  in  terms  of  a  set  of  functions, {xr(r,  £)  }”=i> 
depending  on  nonlinear  parameters,  £  =  (£1,  ...£n, ).  For  a  given  £,  the  optimal  approximation, 
/'( r),  is  the  projection  onto  %: 
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with  the  projector: 
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where  x )  =  ( |  Xi ),  |  X'2  ),  |  Xm  )),  and  the  S  is  the  overlap  matrix  with  elements  Srs  = 

( Xr  |  Xs  )•  The  residual  will  be: 


A2  =  ll  I  f)-&  | /)  ||2  =  </ 1 1  -  ^  | /) 


and  minimizing  it  with  respect  to 

dA2  d 


i  —  1,2,  ...m 


Furthermore: 


i  dX  \  o-i  /  „  i  i  i  „  \  o—i  /  9x  ,  .  ,  .  dS 


-i 


a&  =  la&>s~  (xl  +  lx>ir  (^l  +  lx)  % 


(x 


where 
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(237) 


as 


-i 


d£i 


=  _ s_1  I?  S-1  =  -s 


3-1 


% 


<9x 


<9x 


lx)  +  <x|  ^r) 


di, 


d£i 


3-1 


Replacing  in  eq(237)  and  grouping  terms: 


(238) 


^=l%)S-I(||(l^)  +  (l-^)||)S-I(x 


and  for  real  functions  it  turns: 


(239) 


(240) 


In  general,  A2  will  be  a  complicated  nonlinear  function  of  Eq(240)  can  be  combined  with  a 
gradient  technique  for  minimization,  but  it  is  preferable  to  use  a  second  order  method.  To  do 


this,  the  second  derivatives  of  the  projector  can  be  obtained  from  eq(237 ): 
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Q2g> 

dZfiZi 


d2X 

d£jd£i 


)  s  (x I  +  |x)S 


-1 


d2X  |  |  dx  '  o-i  /  9x 

dtjdti 1  +  1  % ;  '  dti 


and 


dx  \  o-i  /  d*  i  ,  i  d*  \  fdS  /  i,i  v  /^OS  ^  dx 

so)s  <%l  +  l%HArJ<xl  +  lx>UrJ<so 

dx  \  fdS  ^  ,  ,  x  (d S  1  \  ,  dx  , 

so}  l  so  j (xl  +  x>  v  so  J  W 

x>(S)<xl 


<92  A2  <92  & 

%=  aT»T  =  </l  aTarl/) 


dZjdZi 


dZM% 


Notice  that  the  first  derivatives  of  S  1  are  given  by  ec(238 )  and  the  second,  by: 


(241) 

(242) 


d2S 


-i 


9CMt 


__ _d_ 
d£j 

=  S’1 


0—1  0  —  1 

d£i 

dS_  ,  dS  dS  _i  OS  _  02S 

(Kj  0&  +  (Ki  d£j  d£jd£i 


3-1 


where 


d2S  _  /  Q2X  |  v  /  |  d2X  )  (dx,dX)  (dx\dx, 
dZjdZi  v  dZjdZi  lX/^KX  1  dZjdZi 1  K  dZj  '  dZi '  K  0&  1  dZj  1 


As  an  alternative,  eq(  239  >  can  be  used  to  give: 


(243) 


(244) 


* 

+  lx)  I  S-(4%l  (!-<?)  +  (!-<?)  I  Agr}  |  s-1  (X 


dZjdZi 


dZjdZi 


,_i  0%  0^  0^  dx 


x)  1  S  'W  so  so  'so)|Sl<xl 

If  each  approximating  function  depends  on  a  single  parameter  {xs(r,  £s)}£Li : 


(245) 


dXs  _  s  dxi 
dZi  "  dZi 


Defining: 


(246) 


ct  =  (/|x)  =  ((/lxi),(/lx2),...)  c  =  (x|/)  = 


(  X2  |  /  ) 


(247) 
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Da  =  -2  Ci  [C'l  -  (N  C)J  +  [C[  (S~%-  C'  -  67,  T{j  Cj\ 

-  [Cl  (S~%  (MC)j  +  Cj  (S~%  (M  C)j 
+  Cj  (M  S-1)^  Cl  +  Ct  (M  S"%  C'\ 

+  (M  S_1  Cj  +  (MC)i  (S~%-  (M  C)j 

+  C,(MS-%(MC),  +  Q(MS-%(MC),]|  (256) 
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Appendix  D.  Installing  and  running  the  codes 

The  *_global_*  .  f  90  file  containing  the  modules  always  must  be  compiled  first  with  a  suit¬ 
able  FORTRAN  compiler.  In  case  of  the  quadruple  precision  code,  the  compiler  must  admit 
this  possibility. 

Once  the  corresponding  object  and  the  .  mod  files  have  been  generated,  the  file  containing  the 
code  of  the  algorithm  must  be  compiled,  and  the  resulting  objects  must  be  linked. 

In  case  of  the  multiprecision  version,  the  Fortran-90  Multiprecision  System  (mpf  un  90)  must  be 
installed  in  the  computer.  This  package,  created  and  maintained  by  David  H.  Bailey,  is  freely 
distributed  and,  in  the  moment  this  report  is  being  made,  it  can  be  downloaded  in  the  URL: 
http  :  /  / crd.  lbl .  gov/  dhbailey/mpdist  /  Instructions  for  installation  and  usage  of 
the  multiprecision  system  can  be  found  within  the  package. 

The  link  instruction  for  multiprecision  must  include  the  path  to  the  directory  where  the  modules 
specific  of  mpfun90  reside.  Beware  that  these  modules  have  been  created  with  the  same 
compiler  to  be  used  for  compiling  the  sources  for  hybrid  or  exchange  integrals. 

To  run  the  programs,  just  proceed  as  usual:  write  the  executable  name  followed  by  the  inputfile 
and  the  output  file  preceded  by  <  and  >  respectively. 

We  are  also  including  input  samples  and  their  corresponding  output  files  for  testing  that  instal¬ 
lation  has  been  succesful,  and  as  a  guide. 

Hereafter,  some  simple  examples  of  installation  and  run  follow.  In  all  of  them,  we  will  assume 
that  Intel’s  if  ort  FORTRAN90  compiler  will  be  used.  For  other  compilers,  make  the  suitable 
changes  (notice  the  comments  above  for  the  quadruple  precision  version). 

Example  1:  installing  and  running  hybrid_2  010_D 
Type  the  following: 

ifort  -03  -c  hybricL2010_global_D . f 90 

ifort  -03  hybricL2010_D . f 90  hybricL2010_global_D . o  -o  hybricL2010_D 
hybrid_2010_D  <  hybrid_2010_D .  inp  >  outputfile_D 

Compare  the  content  of  output  f  ileJD  with  that  of  hybrid_2  0 10  JD  .  out  coming  with  this 
report. 

Example  2:  instcdling  and  running  hybrid_2010_mp  ( multiprecision ) 

In  this  example,  we  will  suppose  that  mpfun90  has  been  installed  in  the  system  and  that  its 
corresponding  modules  and  objects  reside  in  a  directory  named  /lib/mpfun/f  90/.  Type 
the  following: 

ifort  -03  -c  -I/lib/mpfun90/f  90  hybrid_2 01 0_global_mp . f 90 

ifort  -03  -I/lib/mpfun90/f 90  /lib/mpfun90/f 90/mp* . o  hybrid_2010_global_mp  .  o 
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hybrid_2010_mp  .  f  90  -o  hybrid_2  010_mp 


hybrid_2010 _mp  <  hybrid_2  010_mp .  inp  >  outputf  ile_mp 

Compare  the  content  of  outputf  ile_mp  with  that  of  hybrid_2  010_mp  .  out. 
Further  assistance  can  be  found  in  raf  ael .  lopezguam.  es. 
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Supplementary  material 

A  tar  file  compressed  with  gzip  ( EOARD  093069 Jb.tgz ),  accompanying  this  report,  contains  the 
files  mentioned  in  it. 

Extract  the  content,  with 

tar  -vxzf  EOARD.O  930  6 9_b  .  tgz 

or  alternatively  with 

gunzip  EOARD.O  930  6  9Jd  .  tgz 

followed  by 

tar  -vxf  EOARD.O  930  6 9  Jo  .  tar 

it  will  create  a  directory  named  EOARD  093069  with  the  following  files: 

oneelectrelips-Q.f90 

oneelectrelips-Q.inp 

oneelectrelips-Q.out 

coulomb  32  01 0_shiftop  -D.J90 

coulomb  22  01  Oshiftop  -global  _D.f90 

coulomb 2201 0  shiftop  D.  inp 

coulomb 22010 shiftop  _D.  out 

coulomb 22010 Jntnum2D.f90 

coulomb 220 10-intnum-global  2D. f90 

coulomb 22010 JntnumJD.inp 

coulomb 2201 0  intnumJ).  out 

coulomb 22010 Jransl  2D  ,f90 

coulomb 220 10 Jr  ansi  -global  -D.f90 

coulomb 2201 0 Jransl  JD.  inp 

coulomb 220 10 Jransl -D. out 

coulomb 220 10 Jransl -Q.f90 

coulomb 220 10 Jransl  global  Q.f90 

coulomb 220 10 Jransl -Q.inp 

coulomb 22010 Jransl  -Q.out 

coulomb 2201 0 Jransl  jnp.f90 

coulomb 220 10 Jransl -global  jnp.j90 

hybrid 2201 0_D.f90 

hybrid  2010  global  D.f90 

hybrid 2201 0_D.  inp 

hybrid 2201 0J2>.  out 

hybrid 2201 0-Q.J90 

hybrid 2201 0 -global  _Q.f90 
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hybrid  .201 0_Q.  inp 
hybrid  2201 0_Q.  out 
hybrid 2201 0  mp.J90 
hybrid  J201 0  -global  jnp.f90 
hybrid -201 0_mp.  inp 
hybrid -20 10  jnp.  out 
exchange  2201 0_D.f90 
exchange 2201 0 -globed  -D.f90 
exchange 22 01 0-D.  inp 
exchange 22 01 0_D.out 
exchange 2201 0-Q.f90 
exchange 22 01 0 -globed  -Q.f90 
excheinge 22 01 0_Q.  inp 
excheinge 22010  Q.out 
excheinge 2201 0jnp.f90 
excheinge 22 01 0 -globed  _mp.f90 
excheinge 2201 0  mp.  inp 
excheinge 2201  Ojnp.out 
stngexp22010.f 
stogto22010.f90 
stogto22010.inp 
stogto 22010. out 
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